Using Biot Savart or Ampere's Law you will come to the same problem $B$ is not
defined on the ring.
This is the same problem that trying to find the Electric field $E$ of a puntual charge just in the point where the charge is placed $1/r²$ becomes $\infty$...
You need to use the formula for volumes but using the superficial current density $J$ and integrating on a torus, then the magnetic field is well defined. Notice that:
$\frac{\mu_0}{4\pi}\int_V \frac{Jdv\times r}{|r|^3}=\frac{\mu_0}{4\pi}\int_V \frac{4\pi r² d\Omega dr J\times r}{|r|^3}=\mu_0\int_V \frac{J\times u_r r³ d\Omega dr }{|r|^3}=\mu_0\int_V J\times u_r d\Omega dr$
Thus even if $r \to 0$ the $\infty$ does not appear.
The problem is that solving volume integrals is more complicated that using a line... but in this case I cannot find a better option.
General considerations
I think it will help to carefully write out all the expressions. The PDE should be read component-wise
$$
\partial_{tt}\psi_j(x,t)-\nabla^2\psi_j(x,t)=F_j(x,t)
$$
Where $x\in\mathbb{R}^3$. I've replaced $J$ with a generic source $\mathbf{F}$, $\mathbf{E}$ with a generic field $\boldsymbol{\psi}$, and set all constants to unity. The causal Green's function is
$$
G(r,t)=\frac{\delta(t-r)}{4\pi r}
$$
The solution for $\psi_j$ is
$$
\psi_j(x,t)=\int dt' \int d^3x' \ G(x-x',t-t')F_j(x',t') \ \ + \ \ \text{surface terms}
$$
The surface terms are to match initial conditions which we can specify to vanish (see eg. Zangwill Modern Electrodynamics chapter 20). The delta collapses the $t'$ integral and we are left with
$$
\psi_j(x,t)=\frac{1}{4 \pi}\int d^3x' \ \frac{F_j(x',t-|x-x'|)}{|x-x'|}
$$
Note the arguments of the source appearing in the integrand.
The wave equation for $E$
Because $J$ and $\rho$ form part of a continuity equation, they cannot be specified completely independently. For $\psi=E$, the source should be $F_j=-(\partial_t J_j+\partial_j \rho)$, not just $\partial_t J$. You can check this by deriving it from Maxwell's equations. The expression for $E$ is
$$
E_j(x,t)=-\frac{1}{4\pi} \int d^3x' \ |x-x'|^{-1} \bigg[\partial'_{j} \rho(x',t-|x-x'|) \ + \ \partial'_{t}J_j(x',t-|x-x'|) \bigg]
$$
Note the primes on partial derivatives within the integral: the retarded time is $t_r:=t-|x-x'|$. You can see why it's preferable to write the integral as
$$
E_j(x,t)=-\frac{1}{4\pi} \int d^3x' \ |x-x'|^{-1} \bigg[\partial'_{j} \rho(x',t') \ + \ \partial'_{t}J_j(x',t') \bigg]_{t'=t_r}
$$
With $R_j:=x_j-x_j'$, and using the multivariate chain rule we find (Jackson chapter 6.5)
$$
\big[\partial'_j \rho(x',t') \big]_{t'=t_r} = \partial'_j \big[ \rho(x',t')\big]_{t'=t_r} - \hat{R} \big[ \partial'_t \rho(x',t')\big]_{t'=t_r}
$$
On integrating the term $|x-x'|^{-1}\partial'_j \big[ \rho(x',t')\big]_{t'=t_r}$ by parts, we find Jefimenko's equation for $E_j$.
Finally, if you want to specify a thin wire of arbitrary shape, the expression is not as simple as you may think. Start by parametrizing in $\tau$ such that the wire is the set of points $(x,y,z)=(p_1(\tau),p_2(\tau),p_3(\tau))$, then (up to a proportionality constant)
$$
J_j(\mathbf{x})= \int d\tau \ \frac{dp_j}{d\tau} \delta^{(3)}(\mathbf{x}-\mathbf{p}(\tau))
$$
Best Answer
The magnetic field of a single loop of current is exactly solvable both on- and off-axis, with the horrible disadvantage that the off-axis solution uses elliptic integral functions. The formula you quote seems to be the on-axis field of a single loop.
The exact solution is given in Jackson's Classical Electrodynamics, section 5.5, or in this webpage. From the latter, $$ B_z = \frac{\mu_0 I}{2\pi}\frac{1}{\sqrt{(a+r)^2+z^2}}\left[E(k) \frac{a^2-r^2-z^2}{(a-r)^2+z^2}+K(k)\right] $$ is the on-axis component, while $$ B_\rho = \frac{\mu_0 I}{2\pi}\frac{z/\sqrt{r^2-z^2}}{\sqrt{(a+r)^2+z^2}}\left[E(k) \frac{a^2+r^2+z^2}{(a-r)^2+z^2}-K(k)\right] $$ is the cylindrical polar radial component, where $a$ is the loop radius, $k=\sqrt{\frac{4ar}{(a+r)^2+z^2}}$, and $E(k)$ and $K(k)$ are complete elliptic integrals of the first and second kind, respectively.
This is right horrible but it is suitable for computational purposes. To find the field of a solenoid you should integrate this over an on-axis displacement of the ring. You could attempt this analytically (though I'm not aware of any exact solution) or numerically. The latter makes a bit more sense to me as your solenoid is of course a finite superposition of current loops (though of course those have a finite width).
Added:
The field above is for a single loop in the $x,y$ plane, centred at the origin. To perform the integration, you need to transform this result into the field of a loop around an arbitrary point on the $z$ axis. To do this, you need to change $z$ to the difference $z'=z-b$ and $r$ to the distance $$r'=\sqrt{x^2+y^2+(z-b)^2}=\sqrt{r^2+b^2-2bz}.$$ You also need to transform the field: the component $B_z$ is a cartesian component and transforms well, as does the cylindrical component $B_\rho$, but the spherical radial component $B_r$ from the previous version points away from the centre of the loop and needs to be transformed.
Once this is done, you integrate over $b$ from $-L/2$ to $L/2$. Note that this involves integrating nested square roots inside the elliptical integral functions, so it is unlikely to have a closed-form solution. I would urge you to consider a numerical integration scheme that views your solenoid as a finite (though possibly large) collection of current loops at different displacements $b$ and implementing the sum in your code.
If this is unacceptably complicated, and you have some specific limit in mind, look in Jackson. He gives expressions in the limit $k^2\ll1$ which he claims specialize well to the near-axis ($\theta\ll\pi$), near-the-centre ($r\ll a$) and far-field ($r\gg a$) cases, and do not involve elliptic integrals.
I hope this helps.