CIR model. Is there a closed-form solution or even a good proxy of analytical solution?

Question

Is there a closed-form (analytical) solution for the Cox-Ingersoll-Ross SDE \begin{equation} dr_t=k_r(\theta_r-r_t)dt+\sigma_r\sqrt{r_t}dW_t\tag{1} \end{equation} ?

Notice that $\{r_t\}$ is our process of interest, $k_r$ and $\theta_r$ are constant parameters and $dW_t$ denotes Wiener increment.

I would need a closed-form solution of $(1)$ (or even a good analytical proxy) so as to compare it with solution of my numerical simulation of $(1)$. I am aiming at determining the order of accuracy of my numerical solution, so I would need an analytical solution (or something like that) of $(1)$ so as to compute the average error of my simulation code, that is $\mathbb{E}\left(r_t^{TRUE}-r_t^{SIMULATED}\right)$.

Any suggestion or good source?