Значения для моделирования хаотического рассеяния не соответствуют базовому случаю

Мой первый пост о переполнении стека, будьте осторожны. Я написал код, чтобы отслеживать положение частицы массы M на плоскости x, y на потенциале V (r), описываемом четырехмерной системой уравнений движения.

M(dv/dt)=-grad V(r),      dr/dt=v,

Которые решаются с помощью метода 4-го порядка Рунге-Кутты, где r = (x, y) и v = (vx, vy); теперь состояние частицы определяется x, y и углом theta между вектором v и положительной осью x, где величина скорости определяется выражением

|v|=sqrt(2(E-V(r))/M)

где E - энергия в плоскости, а потенциал V (r) определяется выражением

V(r)=x^2y^2exp[-(x^2+y^2)],

вот код, который я сделал для начальных значений

x(0)=3,
y(0)=0.3905,
vx(0)=0,
vy(0)=-sqrt(2*(E-V(x(0), y(0)))),

где E = 0,260 * (1 / exp (2))

    // RK4
    #include <iostream> 
    #include <cmath> 

    // constant global variables
    const double M = 1.0;
    const double DeltaT = 1.0;

    // function declaration
    double f0(double t, double y0, double y1, double y2, double y3); // derivative of y0
    double f1(double t, double y0, double y1, double y2, double y3); // derivative of y1
    double f2(double t, double y0, double y1, double y2, double y3); // derivative of y2
    double f3(double t, double y0, double y1, double y2, double y3); // derivative of y3
    void rk4(double t, double h, double &y0, double  &y1, double &y2, double &y3); // method of runge kutta 4th order
    double f(double y0, double y1); //function to use

    int main(void) 
    {
      double y0, y1, y2, y3, time, E, Em;
      Em = (1.0/(exp(2.0)));
      E = 0.260*Em;
      y0 = 3.0; //x
      y1 = 0.3905; //y
      y2 = 0.0; //vx
      y3 = -(std::sqrt((2.0*(E-f(3.0, 0.0)))/M)); //vy
      for(time = 0.0; time <= 400.0; time += DeltaT) 
        {
          std::cout << time << "\t\t" << y0 << "\t\t" << y1 << "\t\t" << y2 << "\t\t" << y3 << std::endl;
          rk4(time, DeltaT, y0, y1, y2, y3);
        }


      return 0;
    }

    double f(double y0, double y1)
    {
      return y0*y0*y1*y1*(exp(-(y0*y0)-(y1*y1)));
    }

    double f0(double t, double y0, double y1, double y2, double y3)
    {
      return y2;
    }

    double f1(double t, double y0, double y1, double y2, double y3)
    {
      return y3;
    }


    double f2(double t, double y0, double y1, double y2, double y3)
    {
      return 2*y0*((y0*y0)-1)*(y1*y1)*(exp(-(y0*y0)-(y1*y1)))/M;
    }

    double f3(double t, double y0, double y1, double y2, double y3)
    {
      return 2*(y0*y0)*y1*((y1*y1)-1)*(exp(-(y0*y0)-(y1*y1)))/M;
    }


    void rk4(double t, double h, double &y0, double &y1, double &y2, double &y3) // method of runge kutta 4th order
    {
      double k10, k11, k12, k13, k20, k21, k22, k23, k30, k31, k32, k33, k40, k41, k42, k43;
      k10 = h*f0(t, y0, y1, y2, y3);
      k11 = h*f1(t, y0, y1, y2, y3);
      k12 = h*f2(t, y0, y1, y2, y3);
      k13 = h*f3(t, y0, y1, y2, y3);
      k20 = h*f0(t+h/2, y0 + k10/2, y1 + k11/2, y2 + k12/2, y3 + k13/2);
      k21 = h*f1(t+h/2, y0 + k10/2, y1 + k11/2, y2 + k12/2, y3 + k13/2);
      k22 = h*f2(t+h/2, y0 + k10/2, y1 + k11/2, y2 + k12/2, y3 + k13/2);
      k23 = h*f3(t+h/2, y0 + k10/2, y1 + k11/2, y2 + k12/2, y3 + k13/2);
      k30 = h*f0(t+h/2, y0 + k20/2, y1 + k21/2, y2 + k22/2, y3 + k23/2);
      k31 = h*f1(t+h/2, y0 + k20/2, y1 + k21/2, y2 + k22/2, y3 + k23/2);
      k32 = h*f2(t+h/2, y0 + k20/2, y1 + k21/2, y2 + k22/2, y3 + k23/2);
      k33 = h*f3(t+h/2, y0 + k20/2, y1 + k21/2, y2 + k22/2, y3 + k23/2);
      k40 = h*f0(t + h, y0 + k30, y1 + k31, y2 + k32, y3 + k33);
      k41 = h*f1(t + h, y0 + k30, y1 + k31, y2 + k32, y3 + k33);
      k42 = h*f2(t + h, y0 + k30, y1 + k31, y2 + k32, y3 + k33);
      k43 = h*f3(t + h, y0 + k30, y1 + k31, y2 + k32, y3 + k33);


      y0 = y0 + (1.0/6.0)*(k10 + 2*k20 + 2*k30 + k40);
      y1 = y1 + (1.0/6.0)*(k11 + 2*k21 + 2*k31 + k41);
      y2 = y2 + (1.0/6.0)*(k12 + 2*k22 + 2*k32 + k42);
      y3 = y3 + (1.0/6.0)*(k13 + 2*k23 + 2*k33 + k43);

    }

Проблема здесь в том, что когда я запускаю код с заданными начальными условиями, значения не совпадают с тем, что предполагается, в зависимости от случая, заданного проблемой.

как должен выглядеть рисунок с заданными начальными условиями

Теперь я думаю, что правильно понял реализацию метода, но я не знаю, почему графики не совпадают, потому что, когда я запускаю код, частица уходит от потенциала.

Любая помощь будет оценена по достоинству.