Problems using boundary conditions NDSolve

Question

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?

Answer 1

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