I am using the book "A first course in general relativity" by Bernard Schutz. On page 267 he derives equation 10.54 but leaves out some steps that I am trying to do myself. The following is a picture showing the page.
I tried solving it in MATLAB to avoid calculation error (and also it's faster). My progress so far is:
$$(\rho + p)\frac{d\phi}{dr}=-\frac{dp}{dr}\tag{10.27}$$ Using (10.48) and rearranging gives:
$$\frac{d\phi}{dr}=\frac{4\pi r}{3}\frac{\rho+3p}{1-8\pi r^2\rho /3}.$$
This can now be integrated:
$$ \int^{\phi(R)}_{\phi(0)}d\phi = \int^{R}_{0} \frac{4\pi r}{3}\frac{\rho+3p}{1-8\pi r^2\rho /3} \,dr$$
Where $p$ equals (10.52) note that there is a typo for equation (10.52) it should state "$p$" instead of "$p_c$".
Now, $\phi(R)$ one can easily get from the stated boundary condition $g_{00}(R)=-e^{2\phi(R)}=-(1-2M/R) \Rightarrow \phi(R)=\frac{1}{2}\log(1-\frac{2M}{R})$ (this step I feel confident about being correct since this was in my lecture). However I am a bit unsure on how to achieve $\phi(0)$. A cheaty way is to simply use the result (10.54) which gives (by setting $r=0$) $\phi(0)=\log(\frac{3}{2}(1-\frac{2M}{R})^{1/2}-\frac{1}{2})$. The next issue is that if one performs this integration then one notice that the integral on the right side of the equation gives problem due to solution containing $\log$ which gets negative if one inserts the lower limit $0$. I thought this could be solved by taking the exponent of both sides since the result we want is exp($\phi$) but this also did not work out for me.
I will include by MATLAB code below which I have been using to trying to solve this.
My question is the following. How do I derive the Schwarzschild constant-density interior solution (equation 10.54)?
clc
clear
syms rho M r R phi p
g00 = -(1-2M/R);
p = rho((1-2Mr^2/R^3)^(1/2)-(1-2M/R)^(1/2))/(3(1-2M/R)^(1/2)-(1-2M*r^2 /R^3)^(1/2));
phiR = (1/2)log(1-2M/R);
phi0 = log((3/2)(1-2M/R)^(1/2)-1/2);
LS = 1;
RS = (4pir/3)(rho+3p)(1/(1-8pir^2rho/3)); %10.27
simplify(exp(int(LS,phi,phi0,phiR)))
simplify(exp(int(RS,r)))
