(* sove boundary problem for magnetic *)
(* shielding by spherical shell        *)
(* ----------------------------------- *)

(* phi1 = phi(b<r)    *)
(* phi2 = phi(a<r<b)  *)
(* phi3 = phi(r<a)    *)

phi1=-H0*r*Cos[th]+alpha/r^2*Cos[th]
phi2=beta*r*Cos[th]+gamma/r^2*Cos[th]
phi3=delta*r*Cos[th]

(* bc at outer radius *)
(* ------------------- *)
eq1=(D[phi1,Cos[th]]==D[phi2,Cos[th]])/.r->b

eq2=(mu0*D[phi1,r]==mu*D[phi2,r])/.r->b

(* bc at inner radius *)
(* ------------------- *)

eq3=(D[phi2,Cos[th]]==D[phi3,Cos[th]])/.r->a

eq4=(mu*D[phi2,r]==mu0*D[phi3,r])/.r->a

(* solve bc's *)
(* ----------- *)
sol=Simplify[Solve[{eq1,eq2,eq3,eq4},{alpha,beta,gamma,delta}]]

(* asymptotics as mu->infinity of outer field *)
(* ------------------------------------------- *)
alpha0=alpha/.sol[[1]]

Series[alpha0,{mu,Infinity,0}]

(* asymptotics as mu->infinity for inner field *)
(* -------------------------------------------- *)
delta0=delta/.sol[[1]]

Series[delta0,{mu,Infinity,1}]

(* numerical example *)
(* ------------------ *)
a0=2;
b0=3;
list={H0->1,mu0->1,mu->10,a->a0,b->b0}

alpha1=(alpha/.sol[[1]])/.list
beta1=(beta/.sol[[1]])/.list
gamma1=(gamma/.sol[[1]])/.list
delta1=(delta/.sol[[1]])/.list

phi[x_,y_]:=-1*x+alpha1*x/(x^2+y^2)^(3/2)/; Sqrt[x^2+y^2]>b0
phi[x_,y_]:=beta1*x+gamma1*x/(x^2+y^2)^(3/2)/;(Sqrt[x^2+y^2]>a0 && Sqrt[x^2+y^2]<b0)
phi[x_,y_]:=delta1*x/; Sqrt[x^2+y^2]<a0

M[x_,y_]:={D[phi[x1,y],x1]/.x1->1,D[phi[x,y2],y2]/.y2->y}

g1=StreamPlot[M[x,y],{x,0,4},{y,0,4}]

g2=Plot[{Sqrt[3^2-x^2],Sqrt[2^2-x^2]},{x,0,5},PlotStyle->{{Hue[0.7,0.7,0.7],Thick}},PlotRange->{0,5}]

Show[g1,g2]