I modeled something similar earlier this summer. The method I ended up using is really easy.
Your program needs the inital x, y, and z components and the inital v_x, v_y, v_z components. For each time, you want to find what the acceleration is using a = (GM)/(x+y+z)^2.
Once you have that, you...