#include <iostream>
#include <fstream>
#include <iomanip>
#include <math.h>
#include <windows.h>
using namespace std;


using std::setw;
ofstream out;
double step=0.01;
double n;


double f(double x)
{
    return sqrt(pow(x,2)*sqrt(4*pow(x,2)+9)-3*x*x)/sqrt(2);
}

double df(double x, double y)
{
    return (3*pow(y,3)+6*y*x*x)/(2*x*y*y+3*pow(x,3));
}

double RK_Step(double x, double y)
{
    double f1,f2,f3,f4;
    f1=df(x,y);
    f2=df(x+step/2,y+step/2*f1);
    f3=df(x+step/2,y+step/2*f2);
    f4=df(x+step,y+step*f3);
    y=y+ step/6*(f1+2*f2+2*f3+f4);
    return y;
}

void A_Method(double x, double y)
{
    double num_a[4];
    double dlt=0;
    num_a[0]=df(x,y);//y'
    printf("x \t y \t f(x) \t |y-f(x)|\n");
    for(int i=1; i<4; i++){//
        y=RK_Step(x,y);//
        x+=step;
        num_a[i]=df(x,y);
        out << x << " " << y << " " << f(x) << " " << fabs(y-f(x)) << endl;
        cout << setw(4) << x << " " << setw(8) << y << " " << setw(8) << f(x) << " " << setw(13) << fabs(y-f(x)) << endl;
    }// get 4 y'

    for(int i=3; i<1/step; i++)
    {
        y+=step/24*(55*num_a[3]-59*num_a[2]+37*num_a[1]-9*num_a[0]);
        x+=step;
        for(int j=0; j<3; j++)
        {
            num_a[j]=num_a[j+1];
        }

        num_a[3]=df(x,y);//сдвигаем влево и обновляем последний
        if(dlt<fabs(y-f(x))){
            dlt=fabs(y-f(x));
        }

        out << x << " " << y << " " << f(x) << " " << fabs(y-f(x)) << endl;
        cout << setw(4) << x << " " << setw(8) << y << " " << setw(8) << f(x) << " " << setw(13) << fabs(y-f(x)) << endl;
    }
    out << "maximum delta " << dlt << endl;
    cout << "maximum delta " << dlt << endl;
}


void RK(double x, double y)
{
    double dlt=0;
    printf("x \t y \t f(x) \t |y-f(x)|\n");
    for(int i=0; i<1/step; i++){
        y = RK_Step(x,y);
        x += step;
        out << x << " " << y << " " << f(x) << " " << fabs(y-f(x)) << endl;
        cout << setw(4) << x << " " << setw(8) << y << " " << setw(8) << f(x) << " " << setw(13) << fabs(y-f(x)) << endl;

        if(dlt<fabs(y-f(x))){
             dlt=fabs(y-f(x));
        }
    }
    out << "maximum delta " << dlt << endl;
    cout << "maximum delta " << dlt << endl;
}





int main()
{
    n=1024;
    out.open("Output.txt");
    for(int i=1; i<11; i++){
        step=1/n;
        n/=2;
        out << "Step: " << step << endl;
        cout << "Method R-K:" << endl;
        RK(2,2);
        cout << "*************************************" << endl;
        cout << "A_Method:" << endl;
        A_Method(2,2);
    }
    out.close();
    return 0;
}