  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 scalar field. In each case, long-term stability and second-order
convergence is demonstrated. The case of the massless Klein-Gordon 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.
