Метод Рунге-Кутта 4 порядка для решения системы дифференциальных уравнений

#include <stdio.h>
 
#include <math.h>
#include <iostream>
 
using namespace std;
double f(int i, double x, double y[4]) {
    switch (i) {
    case 1:return y[4]; break;
    case 2: return y[1]; break;
    case 3:return y[2]; break;
    case 4:return y[3]; break;
    default: break;
    }
}
/*************************************************************************
Один шаг метода Рунге-Кутта четвертого порядка для решения
системы дифферециальных уравнений.
 
procedure SystemRungeKuttStep(
const x : Real;
const h : Real;
const n : Integer;
var y : array of Real);
 
Алгоритм совершает один шаг метода для системы
диффуров y[i]'=F(i,x,y) для i=1..n
 
Начальная точка имеет кординаты (x,y[1], ..., y[n])
 
После выполнения алгоритма в переменной y содержится состояние
системы в точке x+h
*************************************************************************/
void step(double x, double h, int n, double y[4])
{
    int i = 0;
    double yt[4];
    double k1[4];
    double k2[4];
    double k3[4];
    double k4[4];
 
    for (i = 1; i <= n; i++)
    {
        k1[i] = h*f(i, x, y);
    }
    for (i = 1; i <= n; i++)
    {
        yt[i] = y[i] + 0.5*k1[i];
    }
    for (i = 1; i <= n; i++)
    {
        k2[i] = h*f(i, x + h*0.5, yt);
    }
    for (i = 1; i <= n; i++)
    {
        yt[i] = y[i] + 0.5*k2[i];
    }
    for (i = 1; i <= n; i++)
    {
        k3[i] = h*f(i, x + h*0.5, yt);
    }
    for (i = 1; i <= n; i++)
    {
        yt[i] = y[i] + k3[i];
    }
    for (i = 1; i <= n; i++)
    {
        k4[i] = h*f(i, x + h, yt);
    }
    for (i = 1; i <= n; i++)
    {
        y[i] = y[i] + (k1[i] + 2.0*k2[i] + 2.0*k3[i] + k4[i]) / 6;
    }
}
/*************************************************************************
Алгоритм решает систему диффуров y[i]'=F(i,x,y) для i=1..n
методом Рунге-Кутта 4 порядка.
 
Начальная точка имеет кординаты (x,y[1], ..., y[n])
 
До конечной точки мы добираемся через n промежуточных
с постоянным шагом h=(x1-x)/m
 
Результат помещается в переменную result[4]
*************************************************************************/
void solvesystemrungekutta(double x, double x1, int steps, double result[4]) {
 
    for (int i = 1; i <= steps - 1; i++)
    {
        step(x + i*(x1 - x) / steps, (x1 - x) / steps, 4, result);
    }
}
 
int main() {
    //первая координата 
    double temporaryresult[4];
    solvesystemrungekutta(2, 3, 3, temporaryresult);
    for (int i = 1; i <= 4; i++) {
        printf("Htpekmnf", temporaryresult[i], '\n');
    }
    char c = getchar();
    return 0;
}