How can I solve iteratively for $\Sigma(r)$, T(r) and h(r)?

Question

My naive attempt:

ClearAll["Global`*"]
Md = 10^(-9); 
P = 10; 
R = 10^4; 
α = 10^(-2); 
ϵ = 10^(-4); 
γ = 10^(-2); 
ke = 0.02*(1 + 0.6625); 
k0 = 5*10^20; 
σ = 5.67*10^-8; 
Rg = 8315; 
c = 3*10^8; 
G = 6.67/10^11; 
M = 2.8*10^30; 

Ωk[r_] := Sqrt[(G*M)/r^3]; 
μ = Md/(3*Pi); 
κ = (27*ke)/(2*σ) Rg/μ; 
Co[r_] := 1; 
β[r_] := 0; 

Do[Σ[
   r_] := κ^(-1/5) α^(-4/5) μ^(3/
      5) Ωk[r]^(2/5)*Co[r]^(-1/5);
 h[r_] := (κ α Σ[r]^2 Ωk[
       r]^(-5) Co[r])^(1/6);
 T[r_] = 1/2 μ/Rg (Ωk[r]^2 h[r]^2)/(1 + β[r]);
 Kkr[r_] := k0 Σ[r]/h[r]*T[r]^(-7/2);
 β[r_] := μ/Rg (4 σ)/(3 c)
    T[r]^3/Σ[r] h[r];

 (*Co[r_] := (1 + β[r])^4*(1 + Kkr[r]/ke)*), {2}]

     Plot[Σ[r],{r,10^4, 10^10}]

     Plot[h[r],{r,10^4, 10^10}]

     Plot[T[r],{r,10^4, 10^10}]

to solve iteratively the system of equations:

$$\begin{array} &\mu=\dfrac{M_d}{3\pi}\\ &\rho(r) = \Sigma(r)\, h(r)\\ &\Omega_K(r)= \sqrt{\dfrac{G M}{r^3}}\\ &\beta(r) = \dfrac{\mu}{Rg}\, \dfrac{4\sigma}{3c}\, \dfrac{T^3(r)}{\rho(r)}\\ &T(r) =\dfrac{1}{2} \, \Omega_k^2(r)\, h^2(r)\, \dfrac{1}{1+\beta(r)}\, \dfrac{\mu}{Rg} \\ &h^6(r)=\kappa \alpha \Sigma^2(r) \, \Omega_k(r)^{-5}C(r)\\ &\Sigma(r) = k^{-1/5}\,\alpha^{-4/5}\mu^{3/5}\, \Omega_k^{2/5}(r)\, C^{-1/5}(r)\\ &C(r)=(1-\beta(r))^4 \left( 1+\dfrac{Kkr}{ke} \right)\\ &Kkr(r)=ko \rho(r) \, T^{-7/2}\\ \end{array}$$

fails. How can I solve it?

The iteration procedure goes like this:

• Impose $C=1$ and $\beta(r)=0$
• Solve for $\Sigma$
• Solve for $T$ and $h$
• Calculate $C$ and $\beta$ and $Kkr(r)$
• Solve for $\Sigma$
• ...

Answer 1

We can use the method of the false transient:

ClearAll["Global`*"]
Md = 10^(-9);
P = 10;
R = 10^4;
α = 10^(-2);
ϵ = 10^(-4);
γ = 10^(-2);
ke = 0.02*(1 + 0.6625);
k0 = 5*10^20;
σ = 5.67*10^-8;
Rg = 8315;
c = 3*10^8;
G = 6.67/10^11;
M = 2.8*10^30;

Ωk[r_] := Sqrt[(G*M)/r^3];
μ = Md/(3*Pi);
κ = (27*ke)/(2*σ) Rg/μ;
ic = {Co[0] == 1, β[0] == 
    0, Σ[
     0] == (μ^(3/5)*Ωk[
         r]^(2/5)) (κ^(-1/5)*α^(-4/5)*Co[0]^(-1/5)),
   h[0] == (κ*α*Σ[0]^2*
       Co[0]/Ωk[r]^5)^(1/6),
   T[0] == (1/2)*(Ωk[r]*h[0])^2*(μ/
       Rg)*(1/(1 + β[0])),
   Kkr[0] == (k0*(Σ[0]*h[0]))/T[0]^(7/2)};

eq = {Σ'[
     t] == -Σ[
       t] + κ^(-1/5) α^(-4/5) μ^(3/
         5) Ωk[r]^(2/5)*Co[t]^(-1/5),
   h'[t] == -h[
       t] + (κ α Σ[t]^2 Ωk[
          r]^(-5) Co[t])^(1/6),
   T'[t] == -T[t] + 
     1/2 μ/Rg (Ωk[r]^2 h[t]^2)/(1 + β[t]),
   Kkr'[t] == -Kkr[t] + k0 Σ[t]/h[t]*T[t]^(-7/2),
   β'[
     t] == -β[t] + μ/
       Rg (4 σ)/(3 c) T[t]^3/Σ[t] h[t],
   Co'[t] == -Co[t] + (1 + β[t])^4*(1 + Kkr[t]/ke)};



t0 = 20; sol = 
 ParametricNDSolveValue[{eq, ic}, {Σ[t], h[t], 
   T[t]}, {t, 0, t0}, {r}]

Table[Plot3D[sol[r][[i]], {t, 0, t0}, {r, 10^4, 10^5}, 
  PlotRange -> All], {i, 3}]

Table[Plot[sol[r][[i]] /. t -> t0, {r, 10^4, 10^5}, 
  PlotRange -> All], {i, 3}]

Answer 2

Not a full answer, just showing how to cast it as a fixed point problem. You have a circular dependency on Co and β, so write this function...

funk[{β_, Co_, r_}] := Module[{Σ, h, T, Kkr, βout, Coout},

  Σ = κ^(-1/5) α^(-4/5) μ^(3/5) Ωk[r]^(2/5)*Co^(-1/5);
  h = (κ α Σ^2 Ωk[r]^(-5) Co)^(1/6);
  T = 1/2 μ/Rg (Ωk[r]^2 h^2)/(1 + β);
  Kkr = k0 Σ/h*T^(-7/2);
  βout = μ/Rg (4 σ)/(3 c) T^3/Σ  h;
  Coout = (1 + βout)^4*(1 + Kkr/ke);

  {βout, Coout, r}
  ];

Now you can run FixedPoint[] or FixedPointList[] and see if it converges. Let it run for 1000 iterations and plot Co. It ends up oscillating between about $10^{54}$ and $0$. The second plot is with the lines through the points.

res = FixedPointList[funk, {.1, 1, 10^4}, 1000];

ListPlot[res[[All, 2]]]
ListLinePlot[res[[All, 2]]]