We present two planet fits to the radial velocity measurements of the M dwarf star GJ 876 that account for the mutual perturbations between the planets, as well as their interactions with the star. We also give the results of long-term integrations testing the stability of some of our best-fit systems. Our Newtonian models fit the data much better than do Keplerian models, reducing the χ2ν from 1.88 to as low as 1.34 for an unstable high-inclination system and 1.43 for a stable low-inclination system. Several different local minima with comparable χ2ν are found in parameter space; thus, our results are not able to provide tight constraints on the inclinations of the orbits to the line of sight and actual planetary masses. Most sets of planetary parameters that we have derived represent systems that are stable for at least 108 yr. Test particles orbiting between the two planets within our best-fit planar system are lost in less than 300 yr.