Is there a way to get values from a set of functions with circular dependencies?

Question

I am using the following code to find iteratively the functions $\Sigma(r)$, $h(r)$ and $T(r)$

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_] := (μ^(3/5)*Ωk[ r]^(2/5))(κ^5^(-1)*α^(4/5)*Co[r]^5^(-1)); 
 h[r_] := (κ*α*Σ[r]* Co[r])/Ωk[r]^5; 
 T[r_] := (1/2)*Ωk[r]* h[r]^2*(μ/Rg)*(1/(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}]

The problem is that the last line Co[r_] := (1 + β[r])^4*(1 + Kkr[r]/ke) makes the kernel crash and I don't understand why.

I am using version 10.0.

Answer 1

I evaluated your code in a clean Mathematica 11.3 notebook. It didn't crash any kernel, but it didn't do any iterative evaluation of Σ[r], h[r] and T[r]. All it did was define those functions and three others two times, giving the same definition each time. That is, the result is same as you would get if you did not wrapped the definitions with Do.

The plots don't work because your definitions of the functions are highly recursive without any code to stop the recursion.

Update

Your functions have the following dependency graph

I really can't see any way to break the complex circularity of the functions. If it is possible, it requires problem domain expertise that I don't (I don't even know what the problems domain is). You will have to use your knowledge of the problem domain or seek the advice of an expert to make your problem tractable.