Abstract We construct a numerical relativity code based on the Baumgarte–Shapiro–Shibata–Nakamura (BSSN) formulation for the gravitational quadratic f ( R ) Starobinsky model. By removing the assumption that the determinant of the conformal 3-metric is unity, we first generalize the BSSN formulation for general f ( R ) gravity theories in the metric formalism to accommodate arbitrary coordinates for the first time. We then describe the implementation of this formalism to the paradigmatic Starobinsky model. We apply the implementation to three scenarios: the Schwarzschild black hole solution, flat space with non-trivial gauge dynamics, and a massless Klein–Gordon (KG) scalar field. In each case, long-term stability and second-order convergence is demonstrated. The case of the massless KG scalar field is used to exercise the additional terms and variables resulting from the f ( R ) contributions. For this model, we show for the first time that additional damped oscillations arise in the subcritical regime as the system approaches a stable configuration.