Spectral Chebyshev Approximation of Cosmic Expansion in $f(R)$ Gravity

We present a numerical framework to study the cosmological background evolution in $f(R)$ gravity by employing a $\mathit{\text{spectral}}\mathit{\text{Chebyshev}}\mathit{\text{collocation}}\mathit{\text{approach}}$. Unlike standard integration methods such as Runge–Kutta that often encounter stiffness and accuracy issues, this formulation expands the normalized Hubble function $E(z)=H(z)\mathrm{/}{H}_0$ as a finite Chebyshev series. The modified Friedmann equation is then enforced at selected Chebyshev–Gauss–Lobatto points, converting the original nonlinear differential equation into a system of algebraic relations for the series coefficients. This transformation yields exponentially convergent and numerically stable solutions over the entire redshift domain, $0<z<z_{max}$, eliminating the need for adaptive step-size control. We apply the method to two widely studied $f(R)$ models, Hu–Sawicki and Starobinsky, and perform a combined analysis using cosmic chronometer $H(z)$ data and the Union~3.0 supernova compilation. The reconstructed expansion histories match observations to within $2\sigma$ over $0<z<2$, producing best-fit parameters of approximately $({\mathrm{\Omega }}_{m0},H_0,{\mathrm{\Lambda }}_{\mathrm{e}\mathrm{f}\mathrm{f}})\simeq (0.29,68,1.2\text{–}2.5H_0^2)$. These results indicates that both models reproduce the observed late-time acceleration while permitting small geometric corrections to $\mathrm{\Lambda }$CDM. Overall, the spectral Chebyshev method provides a precise and computationally efficient framework for probing modified-gravity cosmologies in the precision-data era.