A scale-free phase-field model for martensitic phase transformations (PTs) at finite strains is developed as an essential generalization of small-strain models in Levitas et al. (2004) and Idesman et al. (2005). The theory includes finite elastic and transformational strains and rotations as well as anisotropic and different elastic properties of phases. The gradient energy term is excluded, and the model is applicable for any scale greater than 100 nm. The model tracks finite-width interfaces between austenite and the mixture of martensitic variants only; volume fractions of martensitic variants are the internal variables rather than order parameters. The concept of the effective threshold for the driving force is introduced, which can be either positive (e.g., due to interface friction) or negative (e.g., due to defects and stress concentrators promoting PTs). To reproduce PT conditions obtained from experiments or atomistic simulations under general stress tensor, the effective threshold depends on the stress tensor components. Material parameters are calibrated for martensitic PT between single crystalline cubic Si I and tetragonal Si II phases, which has large transformation strains (εt1 = 0.1753; εt2 = 0.1753; εt3 = −0.447). Finite element algorithms and numerical procedures are implemented in the code deal.II. Multiple 3D problems are solved to study the effect of mesh size, holding time during quasi-static loading, and strain rate on the multivariant microstructure evolution in Si I → Si II PT under uniaxial and hydrostatic loadings. The solution exhibits significant lattice rotations both in Si I and Si II, reproducing the appearance of diffuse grain boundaries in Si I and Si II and transforming them in polycrystals, which corresponds to known experiments. While finer mesh can produce a more detailed microstructure, the solution becomes practically mesh-independent after the mesh size is 80 times smaller than the sample size. When approaching the stationary solution, rough mesh leads to convergence to the correct microstructure faster than the fine mesh, because it neglects fine details in the microstructure. In some regions, reverse PT occurs at continuous compression, despite large transformation hysteresis. For most stationary interfaces, local thermodynamic equilibrium conditions (thermodynamic driving force for the interface motion equal to the effective threshold) are satisfied.