I downloaded your code, but it is too complicate.

It took me a lot of time try to understand it.


The problem I found is :

The evolution equations are:


d(xtest)/dt=vxtest;

d(vxtest)/dt=gy;


However, the gy is defined at constraint page (which is calculated after evolution page was done).

The evolution equations were calculated 4 times (when Runge-Kutta method was selected). -- This has something to do with numerical method.

However, gy was a constant in your code (and it was calculated after evolution page was processed).
You should use user defined function to return gy while evolution was processed.

And all the conditions related gy need to be done inside user's function.

I do not understand what you were doing so I do not know how to help you change your code.