Challenges inherent to high-resolution and high signal-to-noise data as well as model degeneracies can cause systematic biases in analyses of strong lens systems. In the past decade, the number of lens modeling methods has significantly increased, from purely analytical methods, to pixelated and non-parametric ones, or ones based on deep learning. We embraced this diversity by selecting different software packages and use them to blindly model independently simulated Hubble Space Telescope (HST) imaging data. To overcome the difficulties arising from using different codes and conventions, we used the COde-independent Organized LEns STandard (COOLEST) to store, compare, and release all models in a self-consistent and human-readable manner. From an ensemble of six modeling methods, we studied the recovery of the lens potential parameters and properties of the reconstructed source. In particular, we simulated and inferred parameters of an elliptical power-law mass distribution embedded in a shear field for the lens, while each modeling method reconstructs the source differently. We find that, overall, both lens and source properties are recovered reasonably well, but systematic biases arise in all methods. Interestingly, we do not observe that a single method is significantly more accurate than others, and the amount of bias largely depends on the specific lens or source property of interest. By combining posterior distributions from individual methods using equal weights, the maximal systematic biases on lens model parameters inferred from individual models are reduced by a factor of 5.4 on average. We investigated a selection of modeling effects that partly explain the observed biases, such as the cuspy nature of the background source and the accuracy of the point spread function. This work introduces, for the first time, a generic framework to compare and ease the combination of models obtained from different codes and methods, which will be key to retain accuracy in future strong lensing analyses.