I suppose I would need to see the program to answer these questions. I presented my Lagrangian and the solutions based on them. This gives the values in my previous post.
I am not sure where this comes from. The source I have gives a value for GR of:
M(r) = M(∞) / [1 - 3 L^2 / (c^2 r^2)]
Where L is the orbital angular momentum, with L^2 roughly GM(∞)p, where p is the semilatus rectum. (The first order corrective terms for L are not needed here, as this is already a first order correction). This includes the time dilation effects.
http://farside.ph.utexas.edu/teachin...n/node116.html
It looks like you are modeling a Newtonian f this is supposed to be modeling the force, you need to include the increase in energy of Mercury with r. e(r)=e(∞) / [1 - G E(∞) / (2 c^4 r)]^2, so the force is proportional to [1 - G E(∞) / (2 c^4 r)]^-3.
Now, G E(∞)=L^2 c^2/r, and a derivation similar to that in the link below so this really should be M(r)=M(∞) / [1 - L^2 / (2 c^2 r^2)]^3
This should get you half the value.