In the above equation, the second term is a correction due to the fact that
in our model the energy density at the surface of the star,
\(\epsilon _{sup} \equiv \epsilon (P = 0)\)  is not zero [1]}. The functions
\(H\)  and \(\beta \)  can be obtained by solving the following system of differential
equations:
\(H^{\prime }(r) & = & \beta (r), \\\frac{d\beta }{dr} & = & 2\left(1 - \frac{2M(r)}{r} \right)^{-1} H(r)\left\lbrace -2\pi \left[5\epsilon (r) +9P(r)+\frac{\epsilon (r)+P(r)}{dP/d\epsilon }\right]\right. \nonumber \\&&\left. \quad +\frac{3}{r^2} + 2\left(1 - \frac{2M(r)}{r} \right)^{-1}\left( \frac{M(r)}{r^2} +4\pi r P(r) \right)^2\right\rbrace  \nonumber \\&&\quad +\frac{2\beta (r)}{r}\left(1 - \frac{2M(r)}{r} \right)^{-1}\left[-1 + \frac{M(r)}{r} +2\pi r^2\left( \epsilon (r) - P(r) \right)\right].\) 
