Summarizing and expanding on the comments...
First, the given solution only holds in the case $\rho \equiv \rho_0$ is constant (though Wald's wording could stand to be clearer on this point). You won't find a nice general expression for an arbitrary equation of state, which is why numerical results are key for modeling real compact objects like neutron stars.
Proceeding with our assumption, the mass inside a sphere is given by $m = (4\pi/3) \rho_0 r^3$, with the total mass related to the total radius in the same way: $M = (4\pi/3) \rho_0 R^3$. Then the differential equation can be rewritten
$$ \frac{\mathrm{d}P}{\mathrm{d}r} = -(P + \rho_0) \frac{(4\pi/3)\rho_0r^3 + 4\pi r^3P}{r(r - 2(4\pi/3)\rho_0r^3)} = -\frac{4\pi}{3} (P + \rho_0) \frac{(\rho_0 + 3P)r}{1-8\pi\rho_0r^2/3}. $$
We now proceed with separation of variables. The only physics left is in the boundary condition: $P = 0$ at $r = R$.1 With this in mind, we can write
$$ \int_P^0 \frac{\mathrm{d}P}{(P+\rho_0)(3P+\rho_0)} = \int_r^R \frac{r\ \mathrm{d}r}{2\rho_0r^2-3/4\pi}. $$
Both of these integrals are quite solvable using cute tricks from high school calculus, nothing more advanced required.2
The result of integrating is
$$ \frac{1}{2\rho_0} \log\left(\frac{P+\rho_0}{3P+\rho_0}\right) = \frac{1}{4\rho_0} \log\left(\frac{8\pi\rho_0R^2/3-1}{8\pi\rho_0r^2/3-1}\right). $$
Using the earlier definition of $M$ this is clearly equivalent to
$$ \frac{P+\rho_0}{3P+\rho_0} = \left(\frac{1-2M/R}{1-2Mr^2/R^3}\right)^{1/2}, $$
which can be solved for $P$ to yield the desired result.
Finally, it should be noted that ODE solving is not the most interesting thing going on here. What is far more intriguing is how one gets the differential equation to begin with - it is not trivial. Remarkably, Karl Schwarzschild did much of this computation in his 1916 paper "Über das Gravitationsfeld einer Kugel aus inkompressibler Flüssigkeit nach der Einsteinschen Theorie" ("On the Gravitational Field of a Sphere of Incompressible Fluid according to Einstein’s Theory") just a few months after Einstein published his GR field equations.
1 Alternatively one could say $P = P_\mathrm{c}$ at $r = 0$ and integrate the other way, finding a formula for pressure at any radius in terms of the central pressure and mass. The solution can then be inverted to find the radius in terms of the pressure, and the value of this at $P = 0$ would give the size of the object. This procedure make sense if one has a better handle on the microphysics responsible for pressure than on the actual sizes of objects, as is often the case.
2 For the pressure integral, you can complete the square to show that the integrand looks like the derivative of the inverse hyperbolic tangent (or equivalently the inverse tangent with $\sqrt{-1}$ in appropriate places). Then it's a simple matter to invert the exponential formula for $\tanh$ to find the logarithmic version of $\tanh^{-1}$.