This project explores the boundary value problem (BVP) for ordinary differential equations concerning the gravitational field inside a star. The study equates the problem to the electric field inside an atom and reduces the partial differential equation of Poisson’s type to a second order ordinary differential equation, utilising high symmetry. The uniqueness of the solution is ensured by applying two conditions at two ends of the independent variable range. The Numerov’s (Cowell’s) algorithm is employed to solve the equation accurately. However, it is identified that numerical solutions can be very sensitive to the value chosen for the second point, necessitating a recursive scheme. The project also introduces the application of Gauss law, Poisson’s equation, and the Numerov–Cowells algorithm in determining the gravitational potential inside a star given a model radial mass density distribution. The study concludes by discussing the possibility of treating the recursive formula as a tridiagonal system of linear equations and solving it with Gaussian elimination with backward substitution algorithm.