# Problems using boundary conditions NDSolve

NDSolve[
{D[r[dm], dm] == 1/(4π r[dm]^2 ρ2[dm]),
D[P[dm], dm] == (-G dm)/(4π r[dm]^4),
D[L[dm], dm] == ϵ,
P[dm] == (ρ2[dm]/(μ mu ) kb T2[dm])+(a/3 T2[dm]^4),
D[T2[dm], dm] == (3 κ L[dm])/(256 π^2 σ r[dm]^4 T2[dm]^3)},
{r[dm], L[dm], T2[dm], ρ2[dm], P[dm]}, {dm, Mstart, 0.01}]


I want to solve this system of coupled equations, but as you can see it is without boundary conditions. Therefore, when I try to add one like:

L[Mstart] = L0


like so

NDSolve[
{D[r[dm], dm] == 1/(4π r[dm]^2 ρ2[dm]),
D[P[dm], dm] == (-G dm)/(4π r[dm]^4),
D[L[dm], dm] == ϵ,
P[dm] == (ρ2[dm]/(μ mu ) kb T2[dm])+(a/3 T2[dm]^4),
D[T2[dm], dm] == (3 κ L[dm])/(256 π^2 σ r[dm]^4 T2[dm]^3)},
L[Mstart] = L0,
{r[dm], L[dm], T2[dm], ρ2[dm], P[dm]}, {dm, Mstart, 0.01}]


This is still not very well defined, but I get the following error message: (Note that Mstart and L0 are real numbers)

NDSolve::dsfun: L$4740[1.5912*10^30] == 3.846*10^26 cannot be used as a function. The problem with this is that in examples on the reference site they do do it this way, but with me it doesn't seem to work, I'm doing this within a module, so perhaps that has something to do with it? ### Edit The parameters are all numbers, G = 6.674*10^-11; kb = 1.38064852*10^-23; σ = 5.670367*10^-8; mu = 1.66053904*10^-27; μ = 1.28617; ϵ = 400; a=7.56572*10^-16  κ can be set to 1. ### Edit The starting values, or boundary conditions are as follows: L[Mstart]==L0, r[Mstart]==Rstart, \[Rho]2[Mstart]==\[Rho]start, T2[Mstart]==Tstart  With L0=1.5912*10^30 Rstart=5.0112*10^8 \[Rho]=7180.8 Tstart=5.7*10^6  ### Update When I use the method kindly provided by zhk, but with the desired boundary conditions(As given in the second Edit) I recieve the following error: NDSolve::ndsz: At dm$7416 == 1.5879824736744075*^30, step size is effectively zero; singularity or stiff system suspected.

When I added a dynamic [Kappa] it seemed to work a little, but after realizing a mistake, it just gives the same error, so I doubt it has something to do with [Kappa], but for reference, here is the [Kappa] I used:

data1=Import["http://www.uio.no/studier/emner/matnat/astro/AST3310/v18/beskjeder/opacity.txt","Table"];
data2=Drop[data1,{2}];
datalogR=Drop[data2[[1]],1];
datalogT=Drop[data2[[All,1]],1];
datalog\[Kappa]=Drop[Drop[data2,1],0,1];
\[Kappa]function=ListInterpolation[datalog\[Kappa],{{3.75,8},{-8,1}}];

\[Kappa][dm_]=(10^\[Kappa]function[Log[10,T2[dm]],Log[10,\[Rho]2[dm]/(T2[dm]/10^6)^3]])/10


The correction was that I forgot the cube added to the (T2[dm]/10^6) part, perhaps the increased value on [Kappa] was making a possible step size greater then it already is?

• I would be good if you can add the missing parameters, otherwise people will most likely not want to try. – user21 Mar 1 '18 at 14:05
• The problem I see is that you gave the boundary condition as separate argument when it should be included in the 1st argument, the list of equations. – m_goldberg Mar 1 '18 at 15:02

G = 6.674*10^-11;
kb = 1.38064852*10^-23;
\[Sigma] = 5.670367*10^-8;
mu = 1.66053904*10^-27;
\[Mu] = 1.28617;
\[Epsilon] = 400;
Mstart = 0;
L0 = 1;
a = 1;
\[Kappa] = 1;

P[dm_] = (\[Rho]2[dm]/(\[Mu] mu) kb T2[dm]) + (a/3 T2[dm]^4)

Eq1 = D[r[dm], dm] == 1/(4 \[Pi] r[dm]^2 \[Rho]2[dm])

Eq2 = D[P[dm], dm] == (-G dm)/(4 \[Pi] r[dm]^4)

Eq3 = D[L[dm], dm] == \[Epsilon]

Eq4 = D[T2[dm], dm] == (3 \[Kappa] L[dm])/(256 \[Pi]^2 \[Sigma] r[dm]^4 T2[dm]^3)

sol = NDSolve[{Eq1, Eq2, Eq3, Eq4, r[Mstart] == 1, L[Mstart] == L0,
T2[Mstart] == 1, \[Rho]2[Mstart] == 1}, {r[dm], L[dm],
T2[dm], \[Rho]2[dm]}, {dm, Mstart, 0.01}]

Plot[{r[dm], L[dm], T2[dm], \[Rho]2[dm]} /. sol, {dm, Mstart, 0.01}]
`

• Well, at least I know now that it can work, which is a relief, and I forgot to give values for all the "start" values, and the boundary conditions, so I will update my post with them. – Hidde Rinsema Mar 2 '18 at 9:41