Field of a magnetized cylinder (or equivalently, a finite solenoid)

(* exact formula, setting mu0 = 1 and R = 1 *)
b[z_,L_]=1/2 ((z+L/2)/Sqrt[1+(z+L/2)^2]  -  (z-L/2)/Sqrt[1 + (z-L/2)^2])

Plot[b[z,1],{z,-20,20},PlotRange->All]  (* length = 1 x radius *)

Plot[b[z,5],{z,-20,20},PlotRange->All]  (* length = 5 x radius *)

Plot[b[z,20],{z,-20,20},PlotRange->All]  (* length = 20 x radius *)

(* approximate formula, setting mu0 = 1 and R = 1 *)
bapprox[z_,L_]=L/(2 z^3)

(* Some plots of field outside the cylinder, exact and also approximate dipole field *)

Plot[{b[z,1],bapprox[z,1]},{z,0.5,20},PlotRange->{0,.5}]  (* length = 1 x radius *)

Plot[{b[z,5],bapprox[z,5]},{z,2.5,20},PlotRange->{0,.2}]  (* length = 5 x radius *)

Plot[{b[z,20],bapprox[z,20]},{z,10,50},PlotRange->{0,.012}]  (* length = 20 x radius *)