Magnetic Field due to Helmholtz Coils

Clear["Global`*"]
B[x_] := ((μ0 i (R^2))/(2 ((((R^2)+(x^2)))^(3/2)))) (* Field due to a single coil of radius R a distance 'x' away *)

Bleft[x_, s_] := B[x+s/2]
Bright[x_, s_] := B[x-s/2]
Btot[x_, s_] := Bleft[x, s] + Bright[x, s]

Evaluate Derivatives at the center

Simplify[D[Btot[x, R], x]]

% /. x-> 0

Simplify[D[Btot[x, R], {x, 2}]]

% /. x -> 0

Simplify[D[Btot[x, R], {x, 3}]]

% /. x -> 0

Simplify[D[Btot[x, R], {x, 4}]]

% /. x -> 0n(* The 4th derivative is not zero at the center. *)

Plot the Field

Pick some plausible numbers.

R = 0.15; μ0 = 4 π (10^(-7)); i = 5.0;

Btot[x, s]

BHC = Btot[0, R] (* Max value at center of Helmholtz Coil configuration *)

Bmax = 2*B[0] (* Max value for 2 overlapping coils *)

Coil[x_, s_] := Piecewise[{{Bmax/BHC, Abs[Abs[x] - s/2] < 0.01}}, 0] (* Draw coils *)

Manipulate[Plot[{Coil[x*R, s*R], Bleft[x*R, s*R]/BHC, Bright[x*R, s*R]/BHC, Btot[x*R, s*R]/BHC}, {x, -1.5, 1.5}, PlotRange -> {0,Bmax/BHC}, PlotPoints -> 600, ImageSize -> Scaled[0.9], 
LabelStyle->Larger,
AxesLabel->{"x/R", "B / (Subscript[B, HC])"}, PlotLegends -> {"Coils", "Left", "Right", "Total"}], {{s, 2.0}, 0, 2.0, Appearance->"Labeled"}]