Метод Рунге-Кутта для системы ДУ

#include <stdio.h>
#include <iostream>
#include <math.h>
#include <fstream>
using namespace std;
 
const double pi = 3.14159;
 
 double f1(double ksi, double a)  // правая часть моего первого уравнения (энергия)
{
    return((-0.7*atan(3*ksi - 3.8) - 1.3)*sin(a));
}
 
 double f2(double ksi, double a)  // правая часть моего первого уравнения (фаза)
 
{
    return (2 * pi*(1.0 / (0.3*atan(3.1*ksi - 3.8) + 1.3) - a/sqrt(pow(a, 2) - 1));
}
 
int main()
 
{
    int j, N, k;
    double gamma_0, h, pi1, ksi, K1_gamma, K2_gamma, K3_gamma, K4_gamma, K1_phi, K2_phi, K3_phi, K4_phi;
    double Phi[52][80], Gamma[52][80];
 
    gamma_0 = 1.07;
    N = 51;
    pi1 = -pi;
    h = 0.1;
    ksi = 0;
 
    
    
    for (j = 0; j <= N; j++) // задаю параметры для энергии и фазы для всех частиц при кси=0
    {
        Phi[j][0] = pi1 + 2 * pi / 53;
        pi1 = Phi[j][0];
        Gamma[j][0] = gamma_0;
    }
    
    
    for (k = 1; ksi <= 7.8; k++) // сам метод Рунге-Кутта
    {
        for (j = 0; j <= N; j++)
        {
            K1_gamma = f1(ksi, Phi[j][k - 1])*h;
            K1_phi = f2(ksi, Gamma[j][k - 1])*h;
            K2_gamma = f1(ksi + h / 2.0, Phi[j][k - 1] + (h*K1_phi) / 2.0)*h;
            K2_phi = f2(ksi + h / 2.0, Gamma[j][k - 1] + (h*K1_gamma) / 2.0)*h;
            K3_gamma = f1(ksi + h / 2.0, Phi[j][k - 1] + (h*K2_phi) / 2.0)*h;
            K3_phi = f2(ksi + h / 2.0, Gamma[j][k - 1] + (h*K2_gamma) / 2.0)*h;
            K4_gamma = f1(ksi + h, Phi[j][k - 1] + h*K3_phi)*h;
            K4_phi = f2(ksi + h, Gamma[j][k - 1] + h*K3_gamma)*h;
            Gamma[j][k] = Gamma[j][k - 1] + (K1_gamma + 2 * K2_gamma + 2 * K3_gamma + K4_gamma) * 1.0 / 6;
            Phi[j][k] = Phi[j][k - 1] + (K1_phi + 2 * K2_phi + 2 * K3_phi + K4_phi) * 1.0 / 6;
        }
        ksi = ksi + h; // след. шаг итерации
    }
 
    
    for (k = 0; ksi <= 7.8; k++) 
    {
        for (j = 0; j <= N; j++)
        {
            cout << Gamma[j][k] << "           " << ksi << endl;
            cout << Phi[j][k] << "          " << ksi << endl;
        }
    }
 
    system("PAUSE");
}