От конца в начало

void rungeKutta_x(double x, double psi, double u, double t0, double T, double h, std::vector<my_data>& result) {
    double k1[2], k2[2], k3[2], k4[2];
    double dj;
    double t = t0;
    std::vector<std::function<double(double, double, double, double) >> vfunc;
    vfunc.push_back(dxdt);
    vfunc.push_back(dpsidt);
    for (t = t0; t <= T; t += h) {
        for (int j = 0; j < 2; ++j) {
            k1[j] = h * vfunc[j](t, x, psi, u);
        }
 
        for (int j = 0; j < 2; ++j) {
            k2[j] = h * vfunc[j](t + double(0.5)*h, x + double(0.5)*k1[0], psi + double(0.5)*k1[1], u);
        }
 
        for (int j = 0; j < 2; ++j) {
            k3[j] = h * vfunc[j](t + double(0.5)*h, x + double(0.5)*k2[0], psi + double(0.5)*k2[1], u);
        }
 
        for (int j = 0; j < 2; ++j) {
            k4[j] = h * vfunc[j](t + h, x + k3[0], psi + k3[1], u);
        }
 
        x += (k1[0] + double(2.) * k2[0] + double(2.) * k3[0] + k4[0]) / double(6.);
        psi += (k1[1] + double(2.) * k2[1] + double(2.) * k3[1] + k4[1]) / double(6.); 
    }
}