Добавить функционал в готовую программу решения систем дифферинциальных уравнений

double *RungeKutt(int n,double x0,double y0,double h,double (*f)(double x,double y))
{
  register int i;
  double *res,s;
  double k1,k2,k3,k4;
  res=new double [n+1];
  res[0]=y0;
  s=x0;
  for(i=1;i<=n;i++)
  {
    k1=h*f(s,res[i-1]);
    k2=h*f(s+h/2,res[i-1]+k1/2);
    k3=h*f(s+h/2,res[i-1]+k2/2);
    k4=h*f(s+h,res[i-1]+k3);
 
    res[i]=res[i-1]+(k1+2*k2+2*k3+k4)/6;
    s+=h;
  }
  return res;
}

void ProcessRungeKuthMethod3( VectorFn& F, VectorD& Y,
                              double x, double xmax, double h )
{
    VectorD K1, K2, K3, K4; double x0 = x;
    while (x < xmax)
    {
        K1 = h * F( x, Y );
        K2 = h * F( x + h/2, Y + 0.5*K1 );
        K3 = h * F( x + h/2, Y + 0.5*K2 );
        K4 = h * F( x + h, Y + K3 );
        Y += 1./6*( K1 + 2*K2 + 2*K3 + K4 );
        x += h;
        printf("%6.2f    ", x);
        PrintVectorD( Y );
        printf("\n");
    }
}

double **RungeKutt(int n,double x0,double y0,double h,double (*f1)(double x,double y),double (*f2)(double x,double y))
{
  register int i;
  double **res,s;
  double k11,k21,k31,k41;
  double k12,k22,k32,k42;
  res=new double* [n+1];
  for(i=0;i<n+1;i++)
  res[i]=new double[2];
  res[0][0]=x0;
  res[0][1]=y0;
  s=x0;
  for(i=1;i<=n;i++)
  {
    k11=h*f1(s,res[i-1][0]);
    k21=h*f1(s+h/2,res[i-1][0]+k11/2);
    k31=h*f1(s+h/2,res[i-1][0]+k21/2);
    k41=h*f1(s+h,res[i-1][0]+k31);
 
    k12=h*f2(s,res[i-1][1]);
    k22=h*f2(s+h/2,res[i-1][1]+k12/2);
    k32=h*f2(s+h/2,res[i-1][1]+k22/2);
    k42=h*f2(s+h,res[i-1][1]+k32);
 
    res[i][0]=res[i-1][0]+(k11+2*k21+2*k31+k41)/6; // точки по Х
    res[i][1]=res[i-1][1]+(k12+2*k22+2*k32+k42)/6; // точки по У
    s+=h;
  }
  return res;
}

 for(i=0;i<=n;i++)
  {
    Series2->AddXY(t,Yres[i][0],"",Series2->SeriesColor); // Yres[i][0] тут наши решения для х
   Series3->AddXY(t,Yres[i][1],"",Series3->SeriesColor); // Yres[i][1] тут для у
    t+=h;
  }

#include <iostream>
#include <iomanip>
using namespace std;
 
typedef double (*TFunc) (double t, double x, double y);
 
double** RungeKutt(int n, double t0, double x0, double y0, double h, TFunc f1, TFunc f2)
{
    register int i;
    double **res, s;
    double k11, k21, k31, k41;
    double k12, k22, k32, k42;
    res=new double* [n+1];
    for(i=0; i<n+1; i++)
        res[i]=new double[2];
    res[0][0] = x0;
    res[0][1] = y0;
    s = 0;
    for(i=1; i<=n; i++)
    {
        double x = res[i-1][0];
        double y = res[i-1][1];
 
        k11=h*f1(s, x, y);
        k12=h*f2(s, x, y);
 
        k21=h*f1(s+h/2, x+k11/2, y+k12/2);
        k22=h*f2(s+h/2, x+k11/2, y+k12/2);
 
        k31=h*f1(s+h/2, x+k21/2, y+k22/2);
        k32=h*f2(s+h/2, x+k21/2, y+k22/2);
 
        k41=h*f1(s+h,   x+k31,   y+k32);
        k42=h*f2(s+h,   x+k31,   y+k32);
 
        res[i][0] = x + (k11+2*k21+2*k31+k41)/6;
        res[i][1] = y + (k12+2*k22+2*k32+k42)/6;
 
        s+=h;
    }
 
    return res;
}
 
double f1(double t, double x, double y)
{
    return -2*x + 4*y;
}
 
double f2(double t, double x, double y)
{
    return -x + 3*y;
}
 
 
int main()
{
    int n = 60;
    double **res = RungeKutt(n, 0, 0, 1, 0.1, f1, f2);
 
    for (int i=0; i<=n; i++)
    {
        double t = i*0.1;
        cout << setw(3) << t << "  " << setw(12) << res[i][0] << "  " << res[i][1] << endl;
    }
 
}

k32=h*f2(x+k21/2, y+k22/2);