dP/dr = f(r, P(r), M(r))
dM/dr = g(r, P(r), M(r))
P(r0) = P0 > 0
M(x0) = 0
k1 = hf(rn, Pn, Mn)
l1 = hg(rn, Pn, Mn)
k2 = hf(rn + h/2, Pn + k1/2, Mn + l1/2)
l2 = hg(rn + h/2, Pn + k1/2, Mn + l1/2)
k3 = hf(rn + h/2, Pn + k2/2, Mn + l2/2)
l3 = hg(rn + h/2, Pn + k2/2, Mn + l2/2)
k4 = hf(rn + h, Pn + k3, Mn + l3)
l4 = hg(rn + h, Pn + k3, Mn + l3)
k = (1/6) *(k1 + 2k2 + 2k3 + k4)
l = (1/6) *(l1 + 2l2 + 2l3 + l4)
rn+1 = rn + h
Pn+1 = Pn + k
Mn+1 = Mn + l