Solutions for each layer of the shell.  For the purposes of this worksheet, allow space for up to 20 layers.

v=Table[c1[k]*r^(n) + c2[k]*r^(-(n)),{k,1,20}];

dv=D[v,r];

Bn=dv[[1]]/u[1] /. r -> R

Define a function to get B.n just inside the inner boundary of the ABC region; this is called as part of the a minimization that computes permeability.

GetBn[enn_,dee_,seed_]:=Module[{qsub,csub,eqs,m,k,intf,vars},
m=Length[seed];
qsub=Table[u[k]->seed[[k]],{k,1,m}];
qsub=Join[qsub,{R->1,n->enn,d->dee}];
eqs={(v[[1]] /. r->R /. qsub) == 1,(v[[m]] /. r->R+m*d /. qsub) == 0}; (* edges of boundary region *)
For[k=1,k<m,k++, (* interface conditions between layers in the boundary region *)
intf={(dv[[k]]/u[k] /. r -> R+k*d/. qsub) == ( dv[[k+1]]/u[k+1] /. r -> R+k*d/. qsub),
(v[[k]] /. r -> R+k*d/. qsub)==(v[[k+1]] /. r -> R+k*d/. qsub)};
eqs=Join[eqs,intf];
];
vars=Join[Table[c1[k],{k,1,m}],Table[c2[k],{k,1,m}]];
csub=Solve[eqs,vars][[1]];
Bn /. qsub /. csub
]

Now, figure out the permeabilities that satisfy the first m asymptotic boundary conditions.  Do this via an error minimization.  Note that working precision, accuracy, etc., have to be crazy, becuase the calculation isn't well conditioned.

GetPerm[een_,guess_]:=Module[{ee,uu,m,U,seed,Q,k,ic},
m=Length[guess];
ee=Rationalize[een];
seed=Table[U[k],{k,1,m}];
ic=Table[{U[k],guess[[k]]},{k,1,m}];
uu=FindMinimum[Total[Table[(GetBn[k,ee,seed]+k)^2,{k,1,m}]],ic,
WorkingPrecision->1000,MaxIterations->10000,PrecisionGoal->100,AccuracyGoal->100];
N[{uu[[1]],seed /. uu[[2]]}] 
]

Test to see if it works.  Compute the permeabilities for δ = 1 for a 3 rd order BC, assuming all ones as the initial condition.

GetPerm[1,{1,1,1}]

pp={0.001,0.0025,0.005,0.01,0.025,0.05,0.1}

Off[FindMinimum::"precw"];

Off[FindMinimum::"sszero"];

Off[FindMinimum::"fmgz"];

Churn through each order, starting at δ = 1 and reducing δ in small steps until δ = 0.001

MaxOrder=9;
iclist={};
ttt={};
For[nn=1,nn≤MaxOrder,nn++,
(* NewGuess=Table[1,{k,1,nn}]; *)
NewGuess=Join[{1},iclist];
tt={};
For[k=Length[pp],k>0,k--,
NewGuess=GetPerm[pp[[k]],NewGuess][[2]];
tt=Append[tt,{pp[[k]],NewGuess}];
Print[pp[[k]],"  ",NewGuess];
If[k==Length[pp],iclist=NewGuess];
];
ttt=Append[ttt,tt];
]

For[k=1,k≤MaxOrder,k++,Print[MatrixForm[ttt[[k]]]]]

Now, we can create a function that directly gets the collection of permeabilities for any δ and ABC order :

GetPermAlt[dee_,enn_]:=Module[{seed},
seed=Interpolation[ttt[[enn]],InterpolationOrder->1][dee];
GetPerm[dee,seed][[2]]
]

Table of permeabilities needed for 0.1 R stackup case

MatrixForm[Table[{1/(10 k), GetPermAlt[1/(10 k), k]}, {k, 1, MaxOrder}]]