File:Ideal-solenoid-field-3D.png

Bx[x_, z_, R_, L_] = 
  If[Abs[x] < 10^-5, 0, 
   1/(4 Pi L) Sum[(\[Beta] = (x + R)^2 + (z - s L/2)^2;
      m = (4 x*R)/\[Beta];
      s Sqrt[\[Beta]]/
        x ((2 - m) EllipticK[m] - 2 EllipticE[m])), {s, {-1, 1}}]];
Bz[x_, z_, R_, L_] = (x1 = -Abs[x]; 
   2/(4 Pi L) Sum[(\[Beta] = (x1 + R)^2 + (z - s L/2)^2;
      m = 4 x1*R/\[Beta];
      s (z - s L/2)/
        Sqrt[\[Beta]] ((((x1 - R)/(x1 + R) EllipticPi[
              4 x1 R/(x1 + R)^2, m]) /. Indeterminate -> 0) - 
         EllipticK[m])), {s, {-1, 1}}]);
R = 1; L = 4;
plot = Plot3D[
  Sqrt[Bx[x, z, R, L]^2 + Bz[x, z, R, L]^2] L, {x, -3, 3}, {z, -4, 4},
   PlotRange -> {0, 1}, BoxRatios -> {6, 8, 3}, 
  AxesLabel -> {"\[Rho]", "z", "B"}, 
  AxesEdge -> {{-1, -1}, {1, -1}, {1, 1}}, PlotPoints -> 200, 
  MaxRecursion -> 6, ImageSize -> 800, 
  Exclusions -> {{Abs[x] == R, Abs[z] < L/2}}, 
  ExclusionsStyle -> {None, Black}, ClippingStyle -> Black, 
  Mesh -> {11, 15, Range[0, 1, 0.1]}, MeshFunctions -> {#1 &, #2 &, #3 &}, 
  MeshStyle -> {Black, Black, {Gray, Thickness[0.002]}}, 
  ColorFunction -> (Directive[Specularity[0.5, 16], 
      Glow@ColorData["DarkRainbow"][#3]] &), 
  Lighting -> {{"Spot", White, {0, -5, 10}, {Pi/4, 50}}, {"Spot", 
     White, {-9, 5, 4}, {Pi/4, 150}}, {"Spot", 
     White, {-5, -10, 3}, {Pi/4, 100}}, {"Spot", 
     White, {-7, 6, 10}, {Pi/4, 100}}}, BaseStyle -> {FontSize -> 20},
   ViewPoint -> {3, -4, 2}]
SetDirectory[NotebookDirectory[]];
Export["Ideal-solenoid-field-3D.png", plot, ImageResolution -> 300];