I recommend you do the simulations in 2-D, since it's much easier to visualize, and use the following formulas for pairwise potential and force between particles i and j:
potential energy = phi = G*mi*mj*log(rij)where mi and Xi are the mass and 2-D position of particle i, rij=||Xi-Xj|| is the distance between particles i and j, and G is the gravitational constant. Note that the use of the (non-unit) vector Xi-Xj in the numerator makes the force law a 1/r law. This formula is discussed in the paper: Greengard-Rokhlin, J. Comput. Phys. 73, 325-348, 1987, if you're curious. Note that this is not Newton's gravitation law, for which the potential goes as 1/r and the force goes as 1/r^2. This 2-D law tends to give more interesting behavior for 2-D simulations, and it obeys Laplace's equation in 2-D. Whatever potential you choose, make sure that your potential and force are consistent (force = -gradient(phi)).
force = -gradient(phi) = -G*mi*mj*(Xi-Xj)/rij^2
The differential equation says that the mass of particle i (mi) times its acceleration (Xi'') equals the sum of the forces on it, namely the sum of the gravitational forces given by the formula above, due to each of the other particles. Each force is a 2-D vector.
I don't care particularly about units. We can use G=1 and fictitious masses and distances. (Feel free to use more physically accurate quantities if you like, but preserve the qualitative nature of the problem.)
What we have here is a system of ODE's: for n particles it's a system of 2n 2nd order ODE's, or a system of 4n 1st order ODE's. Implement a boundary value solution method -- you can choose from one of the methods discussed in Heath's book. I believe the shooting method will be easiest, since this is a nonlinear differential equation.
Don't proceed until you have the above, and you're convinced (by graphing trajectories) that your moon is in a stable, approximately circular orbit of the chosen radius around the earth, and that the earth is moving only slowly. Note that the coordinate system I've recommended is not fixed relative to the center of mass of the two bodies, so after a long time, both will drift away from the origin.
a) Turn the moon motion into an IVP using the initial velocity from part 1. Then the only BVP is the spaceship. It may require some trial and error to determine the appropriate position of the spaceship at t=tq to put it into orbit. That you can do manually. But solving for the initial velocity of the spaceship, once you've fixed the exact target position of the spaceship at t=tq, should be automated. Note that the moon probably won't pass through the same point at t=tq, because of the moon's gravitational attraction to the spaceship.
b) Keep the moon motion a BVP and add the spaceship BVP, so you'd be solving for both the (new, spaceship-perturbed) moon initial velocity and the spaceship initial velocity.
If you're unable to put your spaceship into orbit around the moon, you can crash it into the moon at t=tq, but there will be a grading penalty, and it will be your job to notify the families of the deceased :-) .
The assignment can be implemented in Matlab or any other language.