Earthquake focal mechanism data provide information about the stress state at the origin of these earthquakes. The inversion methods that are commonly used to infer the stress tensor from focal mechanisms have varying complexity but always rely on a number of assumptions. We present an iterative method built upon a classic linear stress tensor inversion that allows to relax the assumption on shear stress magnitudes while preserving the computational simplicity of the linear problem. Every iteration of our method computes the least-squares solution of the problem, which makes the method fast enough to estimate the inverted parameter errors with non-parametric resampling methods such as bootstrapping. Following previous studies, this method removes the fault plane ambiguity in focal mechanism data by selecting the nodal plane that best satisfies the Mohr-Coulomb failure criterion. We first test the performance and the robustness to noise of the proposed method on synthetic data sets, and then apply it to data from Southern California and from the Geysers geothermal field. We focus the study on investigating the consequences of relaxing the assumption of constant shear stress magnitudes. Our variable shear method successfully generalizes its constant shear counterpart: it is able to perform similarly when the constant shear assumption is a good approximation and provides more accurate results when it is not. We provide the Python package ILSI to implement the proposed method.