Решение системы диффуров методом Рунге-Кутта

#include <stdio.h>
#include <stdafx.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 _tmain(){
    //первая координата 
    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;
}