#include <stdlib.h>
#include <iostream>
#include <conio.h>
using namespace std;


double func1(double y1, double y2){
    return -501*y1+500*y2;
}

double func2(double y1, double y2){
    return 500*y1-501*y2;
}

int main(){
    double y1_0,y2_0,h,tmp1,tmp2;
    double* y1 = new double[4];
    double* y2 = new double[4];
    int it,counter=1;
    cout << "Enter initial conditions: " << endl;
    cout << "y1(0) = " ;
    cin >> y1_0;
    cout << endl << "y2(0) = " ;
    cin >> y2_0;
    cout << endl << "Enter amount of needed itrations: ";
    cin >> it;
    cout << endl << "Enter the step h: ";
    cin >> h;
    y1[0]=y1_0;
    y2[0]=y2_0;
    cout <<"Itration: " <<counter<<"\t" << "y1 = " << y1[0] << "\t" << "y2 = " << y2[0] << endl;
    //Находим первые три точки по методу Эйлера
    for (int i=1; i<3; i++){
        y1[i]=y1[i-1]+h*func1(y1[i-1],y2[i-1]);
        y2[i]=y2[i-1]+h*func2(y1[i-1],y2[i-1]);
        counter++;
        cout <<"Itration: " <<counter<<"\t" << "y1 = " << y1[i] << "\t" << "y2 = " << y2[i] << endl;
    }
    while(counter<it){
        int i=3;
        y1[i]=y1[i-1]+h*( (23.0/12.0)*func1(y1[i-1],y2[i-1]) - (16.0/12.0)*func1(y1[i-1],y2[i-2]) + (5.0/12.0)*func1(y1[i-2],y2[i-2]));
        y2[i]=y2[i-1]+h*( (23.0/12.0)*func2(y1[i-1],y2[i-1]) - (16.0/12.0)*func2(y1[i-1],y2[i-2]) + (5.0/12.0)*func2(y1[i-2],y2[i-2]));
        counter++;
        cout <<"Itration: " <<counter<<"\t" << "y1 = " << y1[i] << "\t" << "y2 = " << y2[i] << endl;

        y1[i-3]=y1[i-2];
        y1[i-2]=y1[i-1];
        y1[i-1]=y1[i];

        y2[i-3]=y2[i-2];
        y2[i-2]=y2[i-1];
        y2[i-1]=y2[i];
    }
    _getch();
    return 0;
}