Накапливается ошибка в решении методом Рунге-Кутты

#include <iostream>
 
using namespace std;
 
double h, a, b, y[4];// y - массив содержащий последовательно x, y и производные y
 
int i, n;
 
double f(double p[4], int j) {
    switch (j)
    {
    case 1:
        return p[2];
    case 2:
        return p[3];
    case 3:
        return 2*p[0]*p[1];
    }
} 
// обратить внимание на порядок уравнений
 
double k1(double p[4], int j) {
    return f(p, j);
}
 
double k2(double p[4], int j) {
    double q[4]; int i;
    q[0] = p[0] + h / 2;
    for (i = 1; i < 4; i++)
    {
        q[i] = p[i] + h / 2 * k1(p, i);
    }
    return f(q, j);
}
 
double k3(double p[4], int j) {
    double q[4]; int i;
    q[0] = p[0] + h;
    for (i = 1; i < 4; i++)
    {
        q[i] = p[i] - h * k1(p, i) + 2 * h * k2(p, i);
    }
    return f(q, j);
}
 
void rungekytt(double p[4]) {
    int j;
    for (j = 1; j  < 4; j ++)
    {
        y[j] = p[j] + h / 6 * (k1(p, j) + 4* k2(p, j) + k3(p, j));
    }
}
 
void main() {
    y[1] = 0.5; y[2] = 1; y[3] = 1;
    a = 0; b = 2; n = 20;
    h = (b - a) / n;
    y[0] = a;
    for (i = 0; i < n; i++)
    {
        cout << y[0] << ' ' << y[1] << endl;
        rungekytt(y);
        y[0] += h;
    }
    cout << y[0] << ' ' << y[1] << endl;
    system("pause");
}