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

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

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

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

}

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

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

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

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