Which coordinate system is most suited for this type of problem? ecliptic? equatorial?
If your simulation is starting from the Earth's surface, or from somewhere close to the Earth (e.g., orbit insertion), you want to use some kind of Earth-centered inertial (ECI) frame initially. At some critical distance from the Moon, your simulation should switch to some kind of Moon-centered inertial (MCI) frame.
Which ECI frame, and which MCI frame, it doesn't matter much. The J2000 frame (very similar to the International Celestial Reference Frame, or ICRF) is widely used. Note that this is not a "true" equatorial frame. It was a "mean" equatorial frame at the J2000 epoch (noon TT on 1 Jan 2000), but it is no longer truly Earth equatorial as the Earth's orientation has precessed and nutated since that time. As the equations of motion take on their simplest form in an inertial (non-rotating) frame, many have migrated to the J2000 frame or the ICRF frame (or the ECI equivalent).
Regarding the critical distance for the switch from ECI to MCI (and if the spacecraft is returning to the Earth, a corresponding switch from MCI to ECI on exceeding a critical distance from the Moon): You will get erroneous results if you don't perform these switches and if you use finite precision arithmetic (e.g.,
double in C, C++ and numpy,
float in python). The switch point has to be somewhere reasonable. Not switching at all is not at all reasonable. Switching from ECI to MCU shortly after translunar injection is not reasonable, nor is switching shortly before lunar landing. The Apollo program used the Moon's gravitational sphere of influence (about 66300 km from the center of the Moon) as the switch point. Another option is the Moon's Hill sphere (about 61600 km from the center of the Moon). Both work.