# Using NSolve and FindRoot

I am attempting to find the coefficients $$c_1,c_2,c_3,c_4,c_5,c_6$$ of the two functions $$u(r) = c_1 r + c_2 (1/r) + c_3I_1(r/l_0) + c_4K_1(r/l_0)$$ and $$\phi(r) = c_5ln(r) + c_6 + (f/a)(\frac{du}{dr} + u(r)/r)$$, where $$I_1$$ and $$K_1$$ are the modified Bessel function of the first and second kind, respectively, $$l_0,f,a$$ are constants. The boundary conditions to determined these coefficients are such that

• $$u(r_i)=u_1$$
• $$u(r_o)=u_2$$
• $$\phi(r_i)=\phi_1$$
• $$\phi(r_o)=\phi_2$$
• $$(\lambda+2\mu)\frac{du}{dr}+\lambda \frac{u}{r}+f(-\frac{d\phi}{dr})$$ at $$r_i$$ = 0
• $$(\lambda+2\mu)\frac{du}{dr}+\lambda \frac{u}{r}+f(-\frac{d\phi}{dr})$$ at $$r_0$$ = 0

The following is my Mathematica code:

Y = 139 10^9;
nu = 0.3;
lambda = (Y*nu)/((1 + nu)*(1 - 2*nu));
mu = 0.5*Y/(1 + nu);
f1 = 1 10^-6; f2 = 1 10^-6; f = f1 + 2*f2; kappa = 1 10^-9;
ls = 0.0 10^-6;
l0 = Sqrt[ls^2 + f^2/((lambda + mu)*kappa)];
dervI1[r_] := 1/2 (BesselI[0, r] + BesselI[2, r]);
dervI0[r_] := 1/2 (BesselI[-1, r] + BesselI[1, r]);
dervI2[r_] := 1/2 (BesselI[1, r] + BesselI[3, r]);
d2I1[r_] := 1/2 (dervI0[r] + dervI2[r]);
dervK1[r_] := -1/2 (BesselK[0, r] + BesselK[2, r]);
dervK0[r_] := -1/2 (BesselK[-1, r] + BesselK[1, r]);
dervK2[r_] := -1/2 (BesselK[1, r] + BesselK[3, r]);
d2K1[r_] := -1/2 (dervK0[r] + dervK2[r]);
u[r_] := c1 r + c2  1/r + c3 BesselI[1, r/l0] + c4 BesselK[1, r/l0];
du[r_] := c1 - c2 1/r^2 + 1/l0 c3 dervI1[r/l0] + 1/l0 c4 dervK1[r/l0];
ddu[r_] := -2 c2 1/r^3 + (1/l0)^2 c3 d2I1[r/l0] + (1/l0)^2 c4 d2K1[r/l0];
\[CurlyPhi][r_] := c5 Log[r] + c6 + f/kappa (du[r] + u[r]/r);
EE[r_] := -(c5 /r) - f/kappa (ddu[r] - u[r]/r^2 + du[r]/r);
\[Sigma]r[r_] := (lambda + 2*mu) du[r] + lambda u[r]/r + (f1 + 2*f2) EE[r];
d\[Sigma][r_] := D[\[Sigma]r[r], r];
eq1 = u[10 10^-6];
eq2 = u[20 10^-6];
eq3 = \[CurlyPhi][10 10^-6];
eq4 = \[CurlyPhi][20 10^-6];
eq5 = l0^2*d\[Sigma][r] /. (r -> 10 10^-6);
eq6 = l0^2*d\[Sigma][r] /. (r -> 20 10^-6);
uin = 0.045 10^-6;
uout = 0.05 10^-6;
\[CurlyPhi]in = 0;
\[CurlyPhi]out = 1;
Sol = NSolve[{eq1 == uin, eq2 == uout, eq3 == \[CurlyPhi]in, eq4 == \[CurlyPhi]out, eq5 == 0, eq6 == 0}, {c1, c2, c3, c4, c5,c6}];
c1 = Sol[[1]][[1]][[2]];
c2 = Sol[[1]][[2]][[2]];
c3 = Sol[[1]][[3]][[2]];
c4 = Sol[[1]][[4]][[2]];
c5 = Sol[[1]][[5]][[2]];
c6 = Sol[[1]][[6]][[2]];
(*root=FindRoot[{eq1==uin,eq2==uout,eq3==\[CurlyPhi]in,eq4==
\[CurlyPhi]out,eq5==0,eq6==0},{{c1,0.0018069515755701755},{c2,2.\
8063598198085684*^-13},{c3,-6.70326128908629*^-14},{c4,-2.\
671175261048626*^-7},{c5,4.016694593136571},{c6,33.86980711923503}},\
AccuracyGoal->Infinity,PrecisionGoal->8]
c1=root[[1]][[2]];
c2=root[[2]][[2]];
c3=root[[3]][[2]];
c4=root[[4]][[2]];
c5=root[[5]][[2]];
c6=root[[6]][[2]];*)
Plot[u[r], {r, 10 10^-6, 20 10^-6}]
Plot[\[CurlyPhi][r], {r, 10 10^-6, 20 10^-6}]


For some other values of $$l_0$$ I can get the correct coefficients, but for the value in the above code, I could not. Should I use different function? I appreciate any help or suggestion. Thank you for reading.

Update

I have solved the problem by using NSolve and SetPrecision

Sol = NSolve[{SetPrecision[eq1 == uin, 500], SetPrecision[eq2 == uout, 500], SetPrecision[eq3 == \[CurlyPhi]in, 500], SetPrecision[eq4 == \[CurlyPhi]out, 500], SetPrecision[eq5 == 0, 500], SetPrecision[eq6 == 0, 500]}, {c1, c2, c3, c4, c5, c6}, Reals];