The presence of scalar fields with non-minimal gravitational interactions of the form \xi |\phi|^2 R ξ | ϕ | 2 R may have important implications for the physics of the early universe. We propose a procedure to solve the dynamics of non-minimally coupled scalar fields directly in the Jordan frame, where the non-minimal couplings are maintained explicitly. Our algorithm can be applied to lattice simulations that include minimally coupled fields and an arbitrary number of non-minimally coupled scalars, with the expansion of the universe sourced by all fields present. This includes situations when the dynamics become fully inhomogeneous, fully non-linear (due to e.g. backreaction or mode rescattering effects), and/or when the expansion of the universe is dominated by non-minimally coupled species. As an example, we study geometric preheating with a non-minimally coupled scalar spectator field when the inflaton oscillates following the end of inflation.