Before answering your question, let me introduce some considerations which are very important for the solution.
While you correctly identify $m$ as the gravitational mass, in General Relativity this is a different concept from the baryonic mass. Roughly speaking, the gravitational mass is the difference between the baryonic one and the binding energy due to gravity.
The second distinction one must be aware of is that $\rho(r)$ appearing in Einstein Equations no longer refers to the baryonic mass density, but encloses the local energy density. In a broad sense, this is a bulk consequence of
$$ E = mc^{2} + \frac{p^{2}}{2m} + \ldots $$
Only in the Newtonian limit we can speak of $\rho$ as the mass density. This situation also have another direct consequence:
- Although $dm/dr = 4\pi r^{2}$ emerges from the mass continuity equation in the Newtonian limit, in General Relativity this equation no longer expresses this concept. As a matter of fact, this equation in GR emerges from
$$G_{00} = \frac{8\pi G}{c^{4}}T_{00}\ .$$
The appropriate general-relativistic version/analogue of mass continuity is given by the conservation of the baryon number 4-vector,
$$\nabla_{\mu}(n_{B}u^{\mu}) = 0\ ,$$
where $n_{B}$ is the baryon number density and $u^{\mu}$ is the fluid 4-velocity. By virtue of "Gauss theorem", for the static and spherically symmetric case this expression can be rearranged in terms of the total number of baryons
$$N_{B} = 4\pi \int^{r}_{0} dx\, x^{2}e^{\Lambda(x)} n_{B}(x)\ ,$$
where we assumed the baryonic number density solely depends on the radial coordinate. We can also turn this into a supplementary differential equation,
$$\frac{dN_{B}}{dr} = 4\pi r^{2}e^{\Lambda(r)}n_{B}(r)\\ $$
$$ = \frac{4\pi r^{2}n_{B}(r)}{\sqrt{1 - \frac{2Gm(r)}{c^{2}r}}}\ .$$
In particular, notice this expression preserves the expected relation $dN_{B} = n_{B}dV$.
- The majority (if not all) Equations of State (EOS) for neutron stars are functional relations of the form $\rho(n_{B})$, $P(n_{B})$. By virtue of the composition of functions and of these being invertible, we can also write $\rho(P)$ and $n_{B}(P)$.
Now, the proper answer. For clarity, let us write down the set of equations describing the structure of a neutron star:
EOS
$$\rho(P(r))\ \ \ ,\ \ n_{B}(P(r))\ .$$
Structure
$$\frac{dm}{dr} = 4\pi r^{2}\rho(r)\\
\frac{dP}{dr} = -\frac{(P(r) + \rho(r) c^{2})e^{2\Lambda(r)}}{r^{2}}\left[\frac{Gm(r)}{c^{2}} + \frac{4\pi r^{3}GP(r)}{c^{4}}\right]$$
Volume and baryon number
$$ \frac{dN_{B}}{dr} = \frac{4\pi r^{2}n_{B}(r)}{\sqrt{1 - \frac{2Gm(r)}{c^{2}r}}}\\
\frac{dV}{dr} = \frac{4\pi r^{2}}{\sqrt{1 - \frac{2Gm(r)}{c^{2}r}}}$$
(Ignore for now the boundary conditions). Treated as functions of $r$, we see the structure equations are independent of $N_{B}(r)$ and $V(r)$ but these later functions explicitly depend of $m(r)$ and implicitly of $P(r)$.
Therefore, you can proceed in two ways for integrating these equations:
Sequentially. First solve the structure, then the volume and the baryon number as now you know $P(r)$ and $m(r)$.
Simultaneously.
In either case, the initial conditions ($r=0$) for $N_{B}$, $V$ and $m$ are the same, i.e. $ = 0$; the central density $\rho_{c}$ determines, throughout the EOS, $P_{c} = P(\rho_{c})$ and $n_{B} = n_{B}(\rho_{c})$. Finally, the radius of the star $r_{\ast}$ is determined as the coordinate where $P(r_{\ast}) = 0$.