This post should have original been the answer to an similar titled question on Mathematica.SE. However, the time I was about to answer, the question has been put on hold by its low apparent quality. So I'll leave it here.

Analysis

Simulation of the magnetic field of a solenoid should be easy. Since a solenoid consists of a series of circular currents. The key of simulation lies in the calculation of magnetic field of a circular current. Then several simple shifts along the axis give all field value generated by individual currents. Adding all values gives the overall field.

Theoretical Work

Circular Current

Based on Biot-Savart law, the magnetic field of a current element \( \mathrm{d}\vec{l} \) is:

$$\mathrm{d}\vec{B}=\frac{\mu0I\mathrm{d}\vec{l}\times\hat{e}r}{4\pi r^2}$$

Since the current has rotational symmetry, so does the field of it. Thus we need only compute the the field on a cross section of the current. Let it be the $$y$$-$$z$$ axis plane.

An analysis of the direction of field shows that, on every point of $$y$$-$$z$$ plane, the field vector is in the plane. Thus for every $$\mathrm{d}\vec{B}$$, only its projection in $$y$$-$$z$$ plane is concerned.

Assume the radius of current $$R$$ and the angular parameter $$\theta$$ starting at $$x$$ axis, increasing by the right-hand rule.

Now compute the field at $$(0,y,z)$$ on the $$x$$-$$z$$ plane.

For an arbitrary current element $$\mathrm{d}l$$ at $$(R,\theta)$$,

r = {0, y, z} - {R Cos[\[Theta]], R Sin[\[Theta]], 0};
dl = {-Sin[\[Theta]], Cos[\[Theta]], 0} R \[DifferentialD]\[Theta];

With Biot-Savart law, we obtain $$\mathrm{d}\vec{B}$$ in component form:

{
  (R z Cos[\[Theta]] \[DifferentialD]\[Theta])/(R^2 + y^2 + z^2 - 2 R y Sin[\[Theta]])^(3/2),
  (R z \[DifferentialD]\[Theta] Sin[\[Theta]])/(R^2 + y^2 + z^2 - 2 R y Sin[\[Theta]])^(3/2),
  (R \[DifferentialD]\[Theta] (R - y Sin[\[Theta]]))/(R^2 + y^2 + z^2 - 2 R y Sin[\[Theta]])^(3/2)
}

The $$x$$ component is irrelevant as discussed above. We can simply drop it in the following numerical computation.

Numerical Computation

Define $$R=1$$.

Let dBfy and dBfz be the differential form of $$y$$ and $$z$$ components of magnetic field, respectively.

dBfy[y_, z_, \[Theta]_] := 
  (R z Sin[\[Theta]])/(R^2 + y^2 + z^2 - 2 R y Sin[\[Theta]])^(3/2);
dBfz[y_, z_, \[Theta]_] := 
  (R (R - y Sin[\[Theta]]))/(R^2 + y^2 + z^2 - 2 R y Sin[\[Theta]])^(3/2);

Let By and Bz be the corresponding integral form of $$y$$ and $$z$$ components of magnetic field.

By[y_, z_] := 
  2 NIntegrate[dBfy[y, z, \[Theta]], {\[Theta], -0.5 \[Pi], 0.5 \[Pi]}];
Bz[y_, z_] := 
  2 NIntegrate[dBfz[y, z, \[Theta]], {\[Theta], -0.5 \[Pi], 0.5 \[Pi]}];

Set a set of sample points:

points = Tuples[Range[-2, 2, 0.2], 2];

The field by a single circular current is easily computed by:

field = Through[{By, Bz}[##]] & @@@ points;

Now we can get some visual impression by:

Thread[{points, field}];
ListStreamPlot[%]

Some warning messages were generated during the plotting, which I think are caused by singularity occurred near the current. But I was too lazy to scrutinize in detail.

The simple idea presented above succeeding, we can simply shift the sample points to generate all necessary field components. I will go for a solenoid with 3 circular currents.

points0 = Tuples[{Range[-2, 2, 0.2], Range[-4, 4, 0.2]}];
points1 = Tuples[{Range[-2, 2, 0.2], Range[-5, 3, 0.2]}];
points2 = Tuples[{Range[-2, 2, 0.2], Range[-3, 5, 0.2]}];
{field0, field1, field2} = (Through[{By, Bz}[##]] & @@@ #) & /@ {points0, points1, points2};
field = field0 + field1 + field2;
data = Thread[{points0, field}];
ListStreamPlot[data]

This is the result:

The complete Mathematica notebook file can be found here.