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$
* ...