Let me try to answer my own question to see if what I am thinking makes sense.
Step 1) There are no free currents, so $\vec{\nabla}\times\vec{H}=0$ and we can define a scalar magnetic potential $U$ so that $\vec{H}=\vec{\nabla}U$. Then it must satisfy $$\nabla^2 U=-\vec{\nabla}\cdot \vec{M}=2ze^{z^2}.$$ Also, $U$ is continuous at the cylinder boundary, and its partial derivatives $\frac{\partial U}{\partial r}$ and $\frac{\partial U}{\partial z}$ are also continuous.
Step 2) Writing the laplacian in cylindrical coordinates and assuming azimuthal symmetry, we have inside the cylinder $$ \frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial U}{\partial r}\right)+\frac{\partial^2 U}{\partial z^2}=2ze^{z^2}.$$ With $U=R_{in}(r)Z_{in}(z)$ we can separate variables to get $$\frac{1}{r}\frac{d}{d r}\left(r\frac{dR_{in}}{d r}\right)=-k^2R_{in}$$ and $$\frac{d^2 Z_{in}}{d z^2}=k^2Z_{in}+2ze^{z^2}Z_{in}$$ for some separation constant $k$.
So the function $R_{in}(r)=J_0(kr)$ is a Bessel function and function $Z_{in}(z)$ is something complicated.
Step 3) Again writing the laplacian in cylindrical coordinates and assuming azimuthal symmetry, we have outside the cylinder $$ \frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial U}{\partial r}\right)+\frac{\partial^2 U}{\partial z^2}=0.$$ With $U=R_{out}(r)Z_{out}(z)$ we can separate variables to get $R_{out}(r)=J_0(qr)$ and now $$\frac{d^2 Z_{out}}{d z^2}=q^2Z_{out},$$ for some separation constant $q$.
Step 4) Continuity of $U$ and its radial derivative imply $q=k$. Continuity of derivative with respect to $z$ for every $z$ is not possible since $Z_{out}$ is very different from $Z_{in}$. I think the only way out is if $U(R,z)=0$ for all $z$, which implies $k$ must be $\gamma_n/R$ where $\gamma_n$ is a zero of $J_0$.
So $$U(r,z)=\sum_n J_0(\gamma_nr/R)Z_{in}(n;z)$$ for $rR$.
Assuming I could find $Z_{in}(n;z)$, I could compute the gradient of this function and it would give me the field $\vec{H}$.
Is this calculation correct? 649