I have written a following mathematica code which involve solution to nonlinear equations. It produce results according to my expectation but it take more than a day in execution which is very annoying. I am not expert in coding but I think my way of solving is inefficient. Can somebody help me what should I change to speedup the execution time without effecting the final results. By the way, I am using mathematica 11.3 on dell core i5, 4 GB ram machine.
1/0 // Quiet;
α = 2.3026*(10^-3);
rho = 2.2*10^-12;
MFD = 4.2;
w0 = MFD/2;
R = List[8, 6, 4, 2, 1];
T0c = 25;
T0 = 25 + 273.15;
surrho = 1.184;
surμ = 1.837*10^-5;
surυ = 1.552*10^-5;
surCp = 1.006*10^3;
surKK = 0.02624;
surg = 9.8;
surβ = 3380*10^-6;
Array[TStack, 5];
For[x = 1, x < 6, x++,
If[R[[x]] > MFD, Rp = MFD, Rp = R[[x]]];
S0 = -(2/(Pi*(2.1*10^-6)^2))*((3.998*10^3)/(Exp[-(2*((62.5/
w0)^2))] - 1));
S = S0*Exp[-2*(r*R[[x]]/w0)^2];
A0 = α*S*10^-18;
surR = R[[x]]*10^-6;
surκ =
surKK/(surCp*
surrho); \
Pr = 0.71;
Tc = T0 + 0.000000002;
TcNew = Tc + 0.000000001;
While[Abs[(((TcNew - Tc)/Tc)*100)] > 10^-10, Print["Tc"];
Print[N[Tc, 25]]; Print["TcNew"]; Print[N[TcNew, 25]]; Tc = TcNew;
Tf = (Tc + T0)/2;
Gr = 8*surg*(surR^3)*(surυ ^-2)*(1/Tf)*(Tc - T0);
Print["Gr"]; Print[Gr];
Nu1 = (0.6 + ((0.387*((Gr*Pr)^(1/
6)))/((1 + ((0.559/Pr)^(9/16)))^(8/27))))^2 ;
Hc = (1/2)*Nu1*surKK*(1/surR); Print["Hc"]; Print[Hc];
ξ = 0.85;
σ = 5.669*10^-8;
Hr = (ξ*σ*((Tc^4) - (T0^4)))/(Tc - T0);
H = Hc + Hr;
H = H*10^-12 ;
KK = (2.14749 - (298.76*Tc^(-1)) + ((20.72*10^3)*
Tc^(-2)) - ((0.54*10^6)*Tc^-3))*10^-6;
Cp = (1/
60.0843)*(55.98 + ((15.40*10^(-13))*Tc) - ((14.4*10^(-5))*
Tc^(-2)));
κ = KK/(rho*Cp);
τ := κ*t;
B0 = A0/KK ;
h = (H/KK);
val = (αn*R[[x]])*BesselJ[1, αn*R[[x]]] - (h*R[[x]])*
BesselJ[0, αn*R[[x]]];
roots = Solve[val == 0 && 0 < αn < 1, αn][[All, 1, 2]];
Print["roots"];
Print[roots];
cn = ((2/(Rp^2))*(1/((BesselJ[
0, (roots[[1 ;; Length[roots]]]*Rp)]^2) + (BesselJ[
1, (roots[[1 ;; Length[roots]]]*Rp)]^2))))^0.5;
b1 = 2/((w0*roots[[1 ;; Length[roots]]])^2);
dn = Abs[
B0*cn*((w0/
2)^2)*(((BesselJ[0, (roots[[1 ;; Length[roots]]]*Rp)]*
Exp[-2*((Rp/w0)^2)]) - 1) +
Sum[(((-1)^l)/((2^(l + 1))*(l!)*Gamma[l + 2]))*((Gamma[l + 1] -
Gamma[l + 1,
b1*((roots[[1 ;; Length[roots]]]*Rp)^2)])/(2*(b1^(l +
1)))), {l, 0, Infinity}])];
Tn = (dn/(roots[[1 ;; Length[roots]]]^2))*(1 -
Exp[-(roots[[1 ;; Length[roots]]]^2)*τ]);
u1 = T0c + (Sum[
Tn[[i]]*cn[[i]]*BesselJ[0, (roots[[i]]*(R[[x]]*r))], {i, 1,
Length[roots]}]) /. t -> Infinity; Print["u1="]; Print[u1];
Print["u1 at r=0"]; Print[N[u1 /. r -> 0]];
Print["u1 at r=R[[x]]"]; Print[N[u1 /. r -> 1]];
TcNew = N[(u1 /. r -> 1 )] + 273.15; Print["Tc="]; Print[Tc];
Print["TcNew="]; Print[TcNew]];
Print["Heat Transfer Co-efficient = "]; Print[H];
TStack[x] = u1; Print["x"]; Print[x];
Print[TStack[x]]]; Print["TStack"]; Print[TStack[1]];
Plot[{TStack[1], TStack[2], TStack[3], TStack[4], TStack[5]}, {r, 0,
1}, Frame -> True, PlotRange -> {{0, 1}, All}, Axes -> {True, True},
BaseStyle -> {FontSize -> 22, FontWeight -> Plain,
FontFamily -> Helvetica},
PlotLegends ->
Placed[LineLegend[(Style[#, 22, Plain,
FontFamily -> Helvetica] &) /@ {"R=8 [μm]",
"R=6 [μm]", "R=4 [μm]", "R=2 [μm]", "R=1 [μm]"},
LegendMarkerSize -> 50, LegendLayout -> {"Row", 3}], {0.7, 0.6}],
AxesOrigin -> {0, 25}, AxesStyle -> Dashed,
FrameLabel -> {"r/R (Normalized Radius)",
"T \!\(\*SuperscriptBox[\([\), \(o\)]\)C]"},
PlotStyle -> {{Orange, Thick}, {Dashed, Red,
Thick}, {Dashing[{0.025, 0.01, 0.025, 0.01}], Purple,
Thick}, {Black, Dashing[{Large}], Thick}, {Blue, Dashing[Tiny],
Thick}}, ImageSize -> 700, Filling -> Axis]
```
