The circular restricted gravitational three-body problem (CRGTBP) is a special case of the gravitational three-body problem which is one of the most important n-body problems. The CRGTBP consists of three bodies; the two bodies of them (which are referred to as the primaries or the primary and the secondary) are moving in a circular orbit about their common mass’s center under the influence of their mutual gravitation and the third body is infinitesimal mass moving under the gravitational significance of the two masses where its mass cannot influence the two masses. Furthermore, the third body has a common plane of movement as defined by both the primary and secondary  - . The problem is tackled by studying the motion of the third body assuming full knowledge of the primary and the secondary with regards to their motions    .
Many physicists and mathematicians have used the method of power series to solve a variety of unsolvable differential equations. The method has been used significantly to solve the celestial mechanic problems as better exactness of inspectional data is required for relatively optimal solutions to the bodies’ dynamical equations of motion. Numerous scientists have used this method including Saad et al.  that used it to get recurrent algorithm for comets under non-gravitational motion. Rabe  used a new computational and iteration method to determine a series of periodic Trojan orbits in the restricted problem of three bodies. Deprit and Price  used the numerical methods to compute characteristic exponents in the planar restricted problem of three bodies. Sharaf et al.  found symbolic solution of the three-dimensional restricted three-body problem and applied it for any given set of initial values.
However, we aim in this paper to employ an approach based on the power series to establish an algorithm or mathematical module using Mathematica to tackle this important problem of circular restricted gravitational dynamical problem. More specifically, we determine the components of the velocity and position vectors of the third body with regard to the CRGTBP.
2. Dynamical Equations of Motion of CRGTBP
The dynamical equations of motion for the CRGTBP are given by the following coupled first-order system as follows 
where , and , are the coordinates of the primary, the secondary and the third body, respectively. More, is the primary’s mass, is the secondary’s mass ( is larger than , where ). is the unit of the force of a gravitational constant; with denoting the distance of the primary and that of the secondary both to the third body     given by
Using Broucke’s method , the following system of first-order differential equations  is thus obtained
3. Solution by Power Series
The equations given in (3) are unsolvable analytically. Therefore, these equations are solved in this section using the power series method and their analytical solutions are given in terms of recurrence relations of the coefficients posed by the variables.
Thus, with the use of the method of power series, the motion’s variables are as follows  
In these power series, the first coefficients are given by the known initial values of as follows 
The following recurrence relations for the remaining coefficients are found using power series 
4. Algorithm for CRGTBP
To generate the components of position and velocity for the third body at any time.
, , , , , , , t and NN.
The components of position and velocity for the third body at any time.
* Module list.
5. Application for CRGTBP
For a numerical example of CRGTBP algorithm, the initial values of components for the position and velocity vectors are considered to be
Then, we get the final value components for the position and velocity vectors as given in Figures 1-12.
It is remarkable to mention here that the accuracy is increased with an increase in the number of terms of the power series, but after NN = 50, the accuracy is fixed. So, we aren’t needed to the more of terms of the power series.
Figure 1. Comparing the exact and approximate values for x when NN = 10.
Figure 2. Comparing the exact and approximate values for y when NN = 10.
Figure 3. Comparing the exact and approximate values for z when NN = 10.
Figure 4. Comparing the exact and approximate values for u when NN = 10.
Figure 5. Comparing the exact and approximate values for v when NN = 10.
Figure 6. Comparing the exact and approximate values for w when NN = 10.
Figure 7. Comparing the exact and approximate values for x when NN = 50.
Figure 8. Comparing the exact and approximate values for y when NN = 50.
Figure 9. Comparing the exact and approximate values for z when NN = 50.
Figure 10. Comparing the exact and approximate values for u when NN = 50.
Figure 11. Comparing the exact and approximate values for v when NN = 50.
Figure 12. Comparing the exact and approximate values for w when NN = 50.
In conclusion, the analytical solutions to the circular restricted gravitational three-body problem (CRGTBP) are determined via the application of the power series method. Also, a module or algorithm is specifically designed and implemented via the help of Mathematica software to find the components of the velocity and position vectors for the third body. Finally, the proposed methodology via the devised module has worked accurately and resulted in reliable results as shown.
 Saad, A.S., Banaszkiewicz, M. and Sitarski, G. (2008) A New Algorithm of the Recurrent Power Series Method for the Non-Gravitational Motion of Comets. Applied Mathematics and Computation, 197, 874-879.
 Deprit, A. and Price, J.F. (1965) The Computation of Characteristic Exponents in the Planar Restricted Problem of Three Bodies, (No. MN-415). Boeing Scientific Research Labs Seattle WA Mathematics Research Lab.
 Sharaf, M.A., Ghoneim, R. and Alshaery, A.A. (2011) Symbolic Solution of the Three Dimensional Restricted Three-Body Problem. Contributions of the Astronomical Observatory Skalnaté Pleso, 41, 1-12.