Численное интегрирование системы дифференциальных уравнений методом Рунге - Кутта

#include "stdafx.h"
#include <iostream>
#include <stdio.h>
#include <conio.h>
#include <math.h>
#include <iomanip>   
#include <fstream>
 
using namespace std;
 
double funk(double y, double* ys, double t)
{
 
    ys[0] = 1;
    ys[1] = 0.1;
    ys[2] = 0;
 
    ys[3] = 5+t*t-2*ys[2]-3*ys[1]-y;
 
    return ys[3];
}
 
double runge_k(double f(double y, double *ys, double t), double *y, double *ys, int n, double tn, double tk, int m, double eps)
{
    double h = (tk - tn)/m;;
    double k1, k2 ,k3, k4;
    double dy;
    
    //m /= 2;
 
    int i = 1;
    double y1 = 0, y2 = 0;
    
    for (int j = 0; j < m; j++)
    {
        k1 = h*f(y[j], ys, tn + j*h);
        k2 = h*f(y[j] + k1/2, ys, tn + j*h/2);
        k3 = h*f(y[j] + k2/2, ys, tn + j*h/2);
        k4 = h*f(y[j] + k3, ys, tn + j*h);
 
        dy = (k1 + 2*k2 + 2*k3 + k4)/6;
        y1 += dy;
    }
 
    
    do 
    {
        if (y2 > 0)
        {
            y1 = y2;
            y2 = 0;
        }
 
        m *= 2;
        h = (tk - tn)/m;
 
        for (int j = 0; j < m; j++)
        {
            k1 = h*f(y[j], ys, tn + j*h);
            k2 = h*f(y[j] + k1/2, ys, tn + j*h/2);
            k3 = h*f(y[j] + k2/2, ys, tn + j*h/2);
            k4 = h*f(y[j] + k3, ys, tn + j*h);
 
            dy = (k1 + 2*k2 + 2*k3 + k4)/6;
            y2 += dy;
        }
 
    }
    while(fabs(y2 - y1/y1) > eps);
 
    return y2;
}
 
int _tmain(int argc, _TCHAR* argv[])
{
    double masy[] = {19.0,18.41,17.84,17.29,16.76,16.25,15.76,15.29,14.84,14.41,14.0,
                    13.61,13.24,12.89,12.56,12.25,11.96,11.69,11.44,11.21,11.0};
 
    double masys[] = {1.0,0.1,0.0,0.0};
 
    double x = 0.0, y = 0.0;
    std::cout << runge_k(funk(y, masys, x), masy, masys, 1.0, 0.0, 2.0, 2.0, 0.0001);
 
    return 0;
}

double runge_k(double (*f)(double , double *, double ), double *y, double *ys, int n, double tn, double tk, int m, double eps)
{}
....
std::cout << runge_k(funk, masy, masys, 1.0, 0.0, 2.0, 2.0, 0.0001);