Schwarzschild constant density interior solution 
I am trying to arrive at the solution $$ e^\phi = \frac{3}{2}(1-2M/R)^{1/2} - \frac{1}{2}(1-2Mr^2/R^3)^{1/2};\quad r \leq R  $$ by solving for $\phi$ in $$ (\rho + p) \frac{d\phi}{dr} = - \frac{dp}{dr};$$  where $\rho $ is constant and provided that $$ \phi(r=R)=\frac{1}{2}\ln(1-2M/R)$$ and $$ p(r) = \rho \frac{(1-2Mr^2/R^3)^{1/2} - (1-2M/R)^{1/2}}{3(1-2M/R)^{1/2} - (1-2Mr^2/R^3)^{1/2}}.$$ 

So I attempted to integrate like this
$$ \int_{0}^{\phi}d\phi = \int_{p(0)}^{p(r)} \frac{-1}{\rho + p} dp$$
$$ e^\phi = \frac{3(1-2M/R)^{1/2} - (1-2Mr^2/R^3)^{1/2}}{3(1-2M/R)^{1/2}-1}.$$
This does not match the correct answer. I believe am missing an integration constant somewhere.
 A: You can not integrate the equations like this. Pressure is not constant and when you tried to change the variables in $$(\rho + P(r))\frac{d \phi}{d r}=-\frac{d P}{d r} \tag{1}$$ form $r$ to $\phi$ you ignored that fact. In $\phi$ it would look like $$(\rho + P(r(\phi))=-\frac{d P}{d \phi},$$
which does not lead anywhere nice because we have no clue how $r(\phi)$ or $P(\phi)$ look like.
What you need to do is integrate (1) in $r\equiv \tilde r$:
 $$\frac{d \phi}{d \tilde r}=-\frac{1}{(\rho + P(\tilde r))}\frac{d P}{d \tilde r} \tag{2}.$$
The RHS of (2) can be integrated when plugging in your explicit expression for $P$ in $r$. When you perform the integral from $r=0$ to $r=\tilde r$ the constant of integration $\phi(0)$ appears on the LHS. This constant can be fixed by the value of $\phi$ on the stellar surface.
This integration in $r$ with the boundary condition at $r=R$ will lead to the desired result.
EDIT: Here the explicit integration: To simplify notation let me introduce the constant $Z\equiv M/R$ and the dimensionless radial variable $s=r/R$. With this variables and OPs relation for the pressure eq. (2) reads 
$$\frac{d \phi}{d s} =\frac{2 s Z}{2 s^2 Z+3 \sqrt{1-2 Z} \sqrt{1-2 s^2 Z}-1}.\tag{3}$$
Integrating the LHS in $s$ is trivial. The RHS is not so easy but with some substitutions, the computer algebra system of onces choosing or Wolfram Alpha it can be integrated to yield:
$$\phi(s) +C = \log[3 \sqrt{1 - 2 Z} - \sqrt{1 - 2 s^2 Z}] \tag{4}.$$
Now imposing the boundary condition for $\phi$ on the stellar surface $s=1$ reveals $C=-\log 2$ and eq. (4) becomes
$$\phi(s) = \log \left[\frac{1}{2} \left(3 \sqrt{1-2 Z}-\sqrt{1-2 s^2 Z}\right)\right] \tag{5},$$
the sought-after result. 
A: By your original differential equation, $\phi = -\ln(\rho+p) + C$.
$e^{\phi} = e^C/(\rho + p)$. Plug in the expression of $p$ and apply the boundary condition, you'll get the solution of $e^\phi$, as simple as that.
