I am very familiar with using Ampere's Law to find the B field inside a solenoid ($\mu_0nI(t)$). Then I can use Faradays to find the E field inside ($\propto\dot I$).
I want to get this result this way using Maxwell's equations inside the solenoid for a harmonic current ($\Re[I_0e^{-i\omega t}]$). Maxwells equations will be $$\nabla\cdot E=\nabla\cdot B=0$$ $$\nabla\times E=i\omega B$$ $$\nabla \times B=-i\omega\mu\epsilon E$$ Which will reduce to $(\nabla^2+k^2)\begin{bmatrix}E\\B\end{bmatrix}=0$. From symmetry, I know that B is only in z and E is only in $\theta$ (cylindrical coords). When I solve the equation for B, since $\nabla\cdot B=0\rightarrow dB/dz=0\rightarrow B(z)=\text{Constant}$
But when I solve helmholtz equation in cylindrical coordinates, I get that $$B=AJ_0(k\rho)$$ SO the field depends radially on location from the center. This does not agree with the constant $\mu nI$. Is this approach correct? If I impose the boundary conditions, (which I am unsure of), will I get something that is not a bessel function?
EDIT: "VEry long" solenoid so edge effects are neglected