How to find the magnetic field of a current using the differential form of Maxwell's equations? To find the magnetic field produced by a long straight wire, one would ise either Biot-Savart law or Ampere's Law in integral form. How do you find this simple result starting from $\nabla \cdot \vec{B} = 0$ and $\nabla \times \vec{B} = \mu_0 \vec{J}$? Let's imagine a current flowing in the $\hat{y}$ direction, then the previous equations are:
$$ \partial_x B_x+ \partial_y B_y + \partial_z B_z = 0$$
$$ \partial_y B_z - \partial_z B_y = 0$$
$$ \partial_x B_y - \partial_y B_x = 0$$
$$ \partial_z B_x - \partial_x B_z = \mu_0 J_y$$
And then? what next?
 A: Edited in the light of first comment below answer
Your problem arises because Maxwell's equations are local equations, applying at a point. Thus the four equations in your question, judging by the last one, apply to a point at which $\mathbf J$ is non-zero (that is a point in your wire). But you wish to find $\mathbf B$ at a point P some distance away from this point. Your local equations, applying to points where $\mathbf J$ is non-zero clearly aren't going to do the job by themselves.
The solution is to apply the local equations to all points between the wire and the field point, P, even though $\mathbf J$, and hence $\nabla \times \mathbf B$, is zero at these intermediate points. You might surround your wire by a net of infinitesimal two-dimensional cells covering the a plane through which the wire passes normally. Imposing the condition that ∇×⃗=0
for all cells except the central, current-carrying one will then give you a line integral around a loop surrounding the wire equal to $\mu_0 I$. I leave the details to you!
The method just described is in fact Stokes's theorem in disguise. I don't think you can do without it. So without apology, here is the slick treatment.
Integrate $$\nabla \times \mathbf B =\mu_0 \mathbf J$$
over a surface S bounded by a closed loop (in your case a circle of radius r centred on the -current-carrying wire), so $$\int_S (\nabla \times \vec B) \cdot d \mathbf S=\mu_0\int_S \mathbf J \cdot d\mathbf S$$
Now apply Stokes's theorem to the left hand side:
$$\int_s \mathbf B \cdot d \mathbf s=\mu_0\int_S \mathbf J \cdot d\mathbf S$$
in which the left hand side is a line integral around the bounding loop.
But we have here Ampère's law, since the right hand side is simply the current, I, through the loop. This is perfectly general, but in your special case, having chosen a loop centred on the current-carrying wire, we know by symmetry that $\mathbf B$ is of equal magnitude, B, all round the loop, so
$$\int_s \mathbf B \cdot d \mathbf s=2 \pi r B$$
So we have
$$B=\frac{\mu_0 I}{2 \pi r}$$
A: As others already pointed out, it is hard
to solve this problem in cartesian coordinates
and from the differential Maxwell equations.
But anyway, here is a rough sketch
without going too much into the details.
The current density $\vec{J}$ is zero everywhere,
except in the wire (at $x=0, z=0$) where it is
infinite and pointing in $\hat{y}$-direction,
in such a way that the total current through
a small circle around the wire is $I$.
This can be described using Dirac delta functions:
$$\vec{J}=\hat{y}I\delta(x)\delta(z)$$
Therefore in your 4th equation you need to write
$I\delta(x)\delta(z)$ instead of $J_y$.
From the symmetry of your situation we try the following:


*

*$B_y$ is zero,

*$B_x$ is independent of $y$,

*$B_z$ is independent of $y$.


Using this approach, from your 4 equations 
the 2nd and 3rd are trivially satisfied.
And the 1st and 4th equation become
$$\begin{align}
\partial_x B_x + \partial_z B_z &= 0 \\
\partial_z B_x - \partial_x B_z &= \mu_0 I\delta(x)\delta(z)
\end{align} \tag{1}$$
By some clever guessing you get the solution
$$\begin{align}
B_x &= -C\frac{z}{x^2+z^2} \\
B_z &= +C\frac{x}{x^2+z^2}
\end{align} \tag{2}$$
with a still unknown pre-factor $C$.
You can easily check the correctness of this solution
by plugging it into the differential equations (1), at
least for outside of the wire ($x\neq 0, z\neq 0$).
For finding the pre-factor $C$ you need to plug the
solution (2) into the second of differential equations (1)
and then integrate it over a small area in the $x$-$z$-plane
containing the wire. Due to the singularity
there (at $x=0, z=0$) this is a tricky business. The
result of this integration is $2\pi C = \mu_0 I$.
So we have the solution
$$\begin{align}
B_x &= -\frac{\mu_0 I}{2\pi}\frac{z}{x^2+z^2} \\
B_z &= +\frac{\mu_0 I}{2\pi}\frac{x}{x^2+z^2}
\end{align} \tag{3}$$
Rewriting this solution (3) from cartesian to cylindrical coordinates gives
$$\vec{B}=\frac{\mu_0 I}{2\pi r}\hat{\phi},$$
where $r$ is the distance from the $y$-axis,
and $\hat{\phi}$ is the azimuthal unit-vector around the $y$-axis.
So finally, we arrived at the same well-known solution
as found by other methods.
A: Is it fine to use Fourier transform?
In general, you can represent the field as
$\vec{B} = \frac{1}{2 \pi} \iiiint\limits_{-\infty}^{+\infty} \vec{ \bf B} e^{i\omega t - i k_x x - i k_y y - i k_z z} dk_x dk_y dk_z d\omega$, 
where $\vec{ \bf B}$ is a Fourier transform for your values
$\vec{ \bf B} =\iiiint\limits_{-\infty}^{+\infty} \vec{B} e^{-i\omega t + i k_x x + i k_y y + i k_z z} dx dy dz dt$.
While $\vec{B}$ is the function of $x, y, z, t$, 
$\vec{ \bf B}$ is just the function of $k_x, k_y, k_z, \omega$.
The same can be done for the current.
Now we substitute this representation into differential equations (ommiting math formalism for changing the order of integration and differentiation), i.e.
$\partial_y B_z - \partial_z B_y = $
$ = \frac{1}{2 \pi} \iiiint\limits_{-\infty}^{+\infty}  \partial_y ({\bf B}_z e^{i\omega t - i k_x x - i k_y y - i k_z z}) dk_x dk_y dk_z d\omega 
- \frac{1}{2 \pi} \iiiint\limits_{-\infty}^{+\infty}  \partial_z ({\bf B}_y e^{i\omega t - i k_x x - i k_y y - i k_z z}) dk_x dk_y dk_z d\omega = $
$ =\frac{1}{2 \pi} \iiiint\limits_{-\infty}^{+\infty} (-ik_y) {\bf B}_z e^{i\omega t - i k_x x - i k_y y - i k_z z} dk_x dk_y dk_z d\omega 
- \frac{1}{2 \pi} \iiiint\limits_{-\infty}^{+\infty}  (-ik_z) {\bf B}_y e^{i\omega t - i k_x x - i k_y y - i k_z z} dk_x dk_y dk_z d\omega = $
$ =\frac{1}{2 \pi} \iiiint\limits_{-\infty}^{+\infty}  (-ik_y {\bf B}_z  + ik_z {\bf B}_y) e^{i\omega t - i k_x x - i k_y y - i k_z z} dk_x dk_y dk_z d\omega$
Well, at the end, in all equations instead of differential operator you get multiplication of integrand on corresponding component of $\vec{k}$. While all the integralas are the same and we must have the uniq solution for the system, we can just equate to zero the integrands instead of initial values, getting
$-ik_x {\bf B}_x - ik_y {\bf B}_y - ik_z {\bf B}_z = 0$
$-ik_y {\bf B}_z + ik_z {\bf B}_y = 0$
$ - ik_x {\bf B}_y + ik_y {\bf B}_x = 0$
$ - ik_z {\bf B}_x + ik_x {\bf B}_z = -\mu_0 ik_y {\bf J}_y$
what is just an algebraic system...
This is valid for any $\vec{J}$. Now you can just put the desired form of the current and find the field. In your case it's, as Thomas wrote, two Delta-functions. Without time dependency we can simplify it to 3-dimensional Fourier. 
A: Biot-Savart law gives us magnetic field completely in steady current situations. So, why do we developed the differential form of equations? The Helmholtz Theorem (it is given in the form at the back of Griffiths Introduction to Electrodynamics) states if we are given the curl and divergence of any vector field $\mathbf F$, i.e $$\nabla \cdot \mathbf F = D $$$$\nabla \times \mathbf F = \mathbf C$$ if $D$ and $\mathbf C$ are known then $\mathbf F$ is known completely. It would become just a matter of solving a partial differential equation (which is a very hard matter to deal with). So, when electrodyamics developed we moved towards the differential representations of equations. Many a times, it is not required (in advanced Electrodyanmics) to know the field explicitly but knowing the divergence and curl confirms (by Helmholtz Theorem) that we know the field and can describe phenomena in electrodynamics.   

To find the magnetic field produced by a long straight wire, one would ise either Biot-Savart law or Ampere's Law in integral form. How do you find this simple result starting from $\nabla \cdot \mathbf B$  and $\nabla \times \mathbf B =\mu_0  \mathbf {J}$ 

We cannot even find magnetic field from Biot-Savart law and Ampere's when situations become complex, they work only in some very symmetrical case. This might help you even more. 
A: It is time to summon the Jefimenko's Equations. Basically these equations are the solutions to the Maxwell's equations. What you need to do is to plug in $J$ and do the integration. See Wikipedia for the equation in its most general form.
In your case $J$ is constant in time and constant in y, and it is just a line current (call it $I$), so the Jefimenko's Equations for magnetic field in your case is reduced to \begin{align}{\mathbf B}((x_0,0,z_0,t)&=\frac{\mu_0}{4\pi}\int_{-\infty}^{\infty} \frac{I(y,t-\sqrt{x_0^2+y^2+z_0^2}/c){\mathbf j}}{({x_0^2+y^2+z_0^2})^\frac{3}{2}}\times (x_0 {\mathbf i}+z_0{\mathbf k})dy \\&=\frac{\mu_0}{4\pi}\int_{-\infty}^{\infty}I\frac{1}{({x_0^2+y^2+z_0^2})^\frac{3}{2}}(-x_0{\mathbf k}+z_0{\mathbf i})dy\\&=\frac{\mu_0}{4\pi}\int_{-\infty}^{\infty}I\frac{\sqrt{x_0^2+z_0^2}}{({x_0^2+y^2+z_0^2})^\frac{3}{2}}\cdot(\mbox{that direction}) dy\\&=\frac{\mu_0 I}{2\pi}\cdot(\mbox{that direction}).\end{align}
It is basically the Biot–Savart law when $J$ is constant. For the integration, I used one of the many on-line integrators. Of course if you manually solve the Maxwell's Equations for your specific case, you end up with Jefimenko's Equations for your specific case. Jefimenko had done that for us so we do not need to do that in the hard way.
