#include "rkf45.c"
#include <iostream>

using namespace std;


int f(int n,double x, double y[],double dydx[])
{
    dydx[0]=-52*y[0]-100*y[1]+exp(-x);
    dydx[1]=y[0]+sin(x);
    return 0;
}

double *mf(double x,double y[])
{
    double *dydx=new double[2];
    dydx[0]=-52*y[0]-100*y[1]+exp(-x);
    dydx[1]=y[0]+sin(x);
    return dydx;
}

int main()
{
    double  h, relerr, abserr, x1, x2;
    int     n=2, flag=1, nfe, maxfe, fail;
    double  y[2], yp[2];

    maxfe  = 5000;
    relerr = 0.0001;
    abserr = 0.0001;
    y[0]=1;
    y[1]=0;
    x1=0;
    x2=2;

    rkfinit(n,&fail);
    cout<< "x1  y[0]            y[1]\n----------------------------------------------------\n";
    for(int step=0; step<20; step++)
    {
        x1=step*0.1;
        x2=x1+0.1;
    rkf45(f,n,y,yp,&x1,x2,&relerr,abserr,&h,&nfe,maxfe,&flag);
    cout<< x1<<"    "<<y[0]<<"      "<<y[1]<<'\n';
    }


    rkfend();

    double zn025[2],
        zn[2],
        zn1[2]={1,0},
        tn,
        tn025,
        t0=0,
        mih=0.001,//в б поменять
        mph=0.1;
    cout<< "\nx1    y[0]            y[1]\n----------------------------------------------------\n";
    for(int step=0; step<(2.0/mph); step++)
    {
        x1=step*mph;
        tn=t0+step*mph;
        tn025=t0+(step+0.25)*mph;

        for(int j=0; j<100;j++)//в б поменять
        {
            zn[0]=zn1[0];
            zn[1]=zn1[1];
            for(int i=0; i<n;i++)
            {
                zn025[i]=zn[i]+mih*0.25*mf(tn,zn)[i];
            }
            for(int i=0; i<n;i++)
            {
                zn1[i]=zn[i]+mih*(-mf(tn,zn)[i]+2*mf(tn025,zn025)[i]);
            }
        }
        cout<< x1+mph<<"    "<<zn[0]<<"     "<<zn[1]<<'\n';
    }





    return 0;
}