>    with(VectorCalculus):
with(plots):

Assume toroidal field can be set up by a radial gradient of a vector potential with only a Z-component

>    SetCoordinates(cylindrical[s,phi,z]):
simplify(Curl(VectorField(<0,0,A(s,z)>)));

Vector(%id = 25524512)

Invert the curl to get a potential that will give the desired toroidal field.

>    assume(r_m<1, r_m>0, R>z, z>0,R>0,R>s);

>    s:='s';

s := 's'

>    B_desired:=s->piecewise(s < r_m*sqrt(R^2-z^2), Bo*s/(r_m), Bo*sqrt(R^2-z^2)*(sqrt(R^2-z^2)-s)/(sqrt(R^2-z^2)-r_m*sqrt(R^2-z^2))):
Az:=s->(-int(B_desired(r), r=0..s)+int(B_desired(r), r=0..sqrt(R^2-z^2)));

Az := proc (s) options operator, arrow; VectorCalculus:-`+`(VectorCalculus:-`*`(VectorCalculus:-int(B_desired(r),r = 0 .. s),-1),VectorCalculus:-int(B_desired(r),r = 0 .. sqrt(VectorCalculus:-`+`(R^2,V...
Az := proc (s) options operator, arrow; VectorCalculus:-`+`(VectorCalculus:-`*`(VectorCalculus:-int(B_desired(r),r = 0 .. s),-1),VectorCalculus:-int(B_desired(r),r = 0 .. sqrt(VectorCalculus:-`+`(R^2,V...

>    display(plot3d(subs(Bo=1,R=1,r_m=.8,Az(s)), s=0..1.1, z=-1..1, grid=[50,50],view=0..0.5, style=patchnogrid),pointplot3d([seq([cos(2*Pi*i/200),sin(2*Pi*i/200),0], i=0..200)], connect=true, color=black));

[Maple Plot]

>    Az(s);

-PIECEWISE([1/2*Bo*s^2/r_m, s <= r_m*(R^2-z^2)^(1/2)],[-1/2*Bo*(2*(R^2-z^2)^(1/2)*s-s^2-r_m*R^2+r_m*z^2)/(-1+r_m), r_m*(R^2-z^2)^(1/2) < s])+1/2*Bo*R^2-1/2*Bo*z^2

>    B_tor:=Curl(VectorField(<0,0,Az(s)>));

B_tor := Vector(%id = 30999220)

>    B_tor2:=simplify(B_tor[1]^2+B_tor[2]^2);

B_tor2 := PIECEWISE([Bo*s/r_m, s <= r_m*(R^2-z^2)^(1/2)],[Bo*(-(R^2-z^2)^(1/2)+s)/(-1+r_m), r_m*(R^2-z^2)^(1/2) < s])^2

>    display(plot3d(subs(Bo=1,R=1,r_m=.8,B_tor[2]), s=0..1, z=-1..1, view=0..1,grid=[50,50],style=patchnogrid),pointplot3d([seq([cos(2*Pi*i/200),sin(2*Pi*i/200),0], i=0..200)], connect=true, color=black));

[Maple Plot]

>    SetCoordinates(polar[r,theta]):
B3D:=MapToBasis(VectorField( <0,B_tor[2]>, 'polar'[s,theta] ), 'cartesian'[x,y]):

>    display(fieldplot(subs({Bo=1,R=1,r_m=.8,z=.6},B3D), x=-1..1, y=-1..1), plot([cos(theta), sin(theta), theta=0..2*Pi]));
plot([subs({x=0,Bo=1,R=1,r_m=.8,z=0},B3D[1]^2)], y=0..1);

[Maple Plot]

[Maple Plot]

Now for the poloidal field, we assume we can use a toroidal vector potential.  Since we want a dipole type field, we set A_phi=A(r)*sin(theta)

>    SetCoordinates(spherical[r,theta,phi]):
simplify(Curl(VectorField(<0,0,Ap(r)*sin(theta)>)));

Vector(%id = 31265660)

Since we need A_phi to go to zero at r=0 and we need A_phi to go to zero as well as B_pol to go to zero at r=R-r , we try

>    Ap:=Bo/(2*R^2)*r*(R-r)^2;

Ap := 1/2*Bo/R^2*r*(R-r)^2

>    plot(subs(Bo=1,R=1,Ap(r)), r=0..1);

[Maple Plot]

>    B_pol:=simplify(Curl(VectorField(<0,0,Ap*sin(theta)>)));

B_pol := Vector(%id = 31265700)

>    B_pol2:=simplify(B_pol[1]^2+B_pol[2]^2);

B_pol2 := Bo^2*(R-r)^2*(2*cos(theta)^2*R*r-3*cos(theta)^2*r^2+R^2-4*R*r+4*r^2)/R^4

>    SetCoordinates(polar[r,theta]):
B3D:=MapToBasis(VectorField( <B_pol[1],B_pol[2]>, 'polar'[r,theta] ), 'cartesian'[x,y]):
display(fieldplot(subs({Bo=1,R=1},B3D), x=-1..1, y=-1..1), plot([cos(theta), sin(theta), theta=0..2*Pi]));
plot([subs({x=0,Bo=1,R=1},B3D[1]^2)], y=0..1);

[Maple Plot]

[Maple Plot]

>