I believe you have a mistake in your formula as the self-inductance of a coil is given by $$L\approx\mu_0 \frac{n^2 A}{\ell};$$
here $n$ is the number of windings, $A$ is area of the cross-section, and $\ell$ is the length of the coil.
Your task is to maximize $L$ with the constraint that the length of the copper wire is $W$. Assuming that the solenoid is a cylinder, the cross-section read $A=\pi R^2$ with $R$ the radius of the cylinder.
A solenoid with $n$ windings needs a wire of length $W= 2\pi Rn$. Thus,
$$ L \approx \mu_0 \frac{W^2}{\ell}.$$
We see that the inductance of the solenoid decreases with increasing length (keeping the total length of the wire fixed). Thus, we obtain the largest self-inductance having the smallest length which is a single loop with $n=1$. For a single loop the formula given above is not correct (as it assumes $\ell \gg \sqrt{A}$) and thus we have
$$L\approx \mu_0 R \ln (R/r) \approx \mu_0 \frac{W}{2\pi} \ln(W/r) $$
with $r$ the radius of the wire.
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.
Best Answer
$n$ in this case isn't the number of loops, it's the number of loops per unit length. This is why you have to multiply it by $L$ to get $N$.