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

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$$
• ...
• Does FindRoot not find the solution when you give it decent starting values? – Sjoerd Smit Jan 31 at 14:31
• @SjoerdSmit I need to plot the functions $\Sigma$ $h$ and $T$. – mattiav27 Jan 31 at 14:40
• It should possible to use FindRoot for a range of values of r, solving the equations for each individual value, right? In fact, once you have a solution for 1 value of r, you could use that solution as an initial guess for neighbouring values of r. – Sjoerd Smit Jan 31 at 15:09
• Code and equations do not match. – Alex Trounev Jan 31 at 15:09
• @AlexTrounev I ìt should be fine now – mattiav27 Jan 31 at 15:13

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 == 1, β ==
0, Σ[
0] == (μ^(3/5)*Ωk[
r]^(2/5)) (κ^(-1/5)*α^(-4/5)*Co^(-1/5)),
h == (κ*α*Σ^2*
Co/Ωk[r]^5)^(1/6),
T == (1/2)*(Ωk[r]*h)^2*(μ/
Rg)*(1/(1 + β)),
Kkr == (k0*(Σ*h))/T^(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}] • What does t do'? If I change it, the shape of the functions change. – mattiav27 Aug 21 at 7:04
• @mattiav27 t is an iterative parameter. Here we use a continuous parameter in the model, but when numerically integrating a system of equations, a discrete parameter is automatically used. – Alex Trounev Aug 21 at 9:51

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]]]  