rmax = 100, g = 3, n = 10^10, B = 160^4, T = 170, r0 = 1.477, L=n/(4/3*pi*r^3)
a = 1/(20*r0^3), mu = ((6*L*pi^2)/g)^1/3,
Er® = 3*pr(r)+4B-16/27*mu*T(mu^2-2*mu*T+T^2)

s =
  NDSolve[
    {m'[r] == a*r^2 Er[r],
     pr'[r] == -(r0/r)*((pr[r] +Er[r]) (m[r] + a*r^3* pr[r]))/(r - 2 r0*m[r]),
     m[1] == 1, pr[1] == 1},
    {m, pr}, {r, 1, rmax}]