I am dealing with some model describing the given composition of a white dwarf made with photons, degenerate electrons, positrons and ions. From that, I can get tabulated values of densities ranging from small to large numbers and corresponding to it the total pressure and total energy. Now I'm interested in getting the parameters like density, mass and pressure with respect to the radius at some given initial conditions. I can do it using two coupled ordinary differential equations, i.e.: $$ \frac{\partial m}{\partial r}=4\pi r^2\rho(r)\\ \frac{\partial P}{\partial r}=-\frac{Gm(r)\rho(r)}{r^2} $$ This can be solved, for example, using the forward Euler method which will give us as follows: $$ m(i)=m(i-1)+4\pi r^2 \rho(i-1)\Delta r\\ P(i)=P(i-1)-\frac{Gm(i-1)\rho(i-1)}{r^2}\Delta r $$ Firstly, I want to implement the following initial conditions i.e. $$ m(0)=0 \ \ ; \ \ \rho(0)=\rho_{central} $$ I'm going outward with my calculations. Using my tabulated values of pressure and density I can easily calculate for small $\Delta r$ how the mass and pressure change with respect to the radius breaking the computations when the density reaches a very small value (the surface of the WD). Having all of these I can determine what is the mass of the WD for a given central density.
Secondly, I want to use different initial conditions to get, in the end, the central density for a given mass, i.e. $$ m(R)=M \ \ ; \ \ \rho(R)=0 \ \ or/and \ \ P(R)=0 $$ It means that I have to do my calculations inward. In that case, I have to modify the above result of the Euler method to get: $$ m(i-1)=m(i)-4\pi r^2 \rho(i)\Delta r\\ P(i-1)=P(i)+\frac{Gm(i)\rho(i)}{r^2}\Delta r $$ And this is the moment I got stuck for several weeks. If I take my initial conditions which are the total mass, density and pressure are equal to zero at the surface then the only parameter that stops me from proceeding with these computations is the total radius. And here I need some help or hint. How can I perform these calculations inward when I don't know this radius? Is there maybe some way of transformation of getting rid of this unknown? Or perhaps I can do computations outward even though I've got surface boundary conditions?