#include <stdio.h>
#include <math.h>
#include <conio.h>
#include <iostream>
using namespace std;
#define E 1e-6

void Eiler(double, double, double, double);
void Runge_Kutta_4(double, double, double, double);
void AdamsMeth(double a, double b, double y0, double h, double exp);

int main()
{
	double a = 0.7, b = 1.7, h = 1e-1, y0 = 2.1;
	Eiler(a, b, h, y0);
	cout << endl;
	Runge_Kutta_4(a, b, h, y0);
	AdamsMeth(a, b, y0, h, E);
	_getch();
	return 0;
}

double function(double X, double Y)
{
	return X + cos(Y / 1.25);
}

void Eiler(double A, double B, double h, double y)
{
	double X[11], Y[11];
	Y[0] = y;
	for (int i = 0; i < 11; i++)
	{
		X[i] = A + i * h;
	}
	printf_s("Eiler method:\ni = %d, X = %lf, Y = %lf\n", 0, X[0], Y[0]);
	for (int i = 1; i < 11; i++)
	{
		Y[i] = Y[i - 1] + function(X[i - 1], Y[i - 1])*h;
		printf_s("i = %d, X = %lf, Y = %lf\n", i, X[i], Y[i]);
	}
	//printf_s("Answer: %lf\n", function(X[10], Y[10]));
}

void Runge_Kutta_4(double A, double B, double h, double y)
{
	double X[11], Y[11], dY[10], k1, k2, k3, k4;
	Y[0] = y;
	for (int i = 0; i < 11; i++)
	{
		X[i] = A + i * h;
	}
	printf_s("Runge-Kutta method:\ni = %d, X = %lf, Y = %lf\n", 0, X[0], Y[0]);
	for (int i = 1; i < 11; i++)
	{
		k1 = h * function(X[i - 1], Y[i - 1]);
		k2 = h * function(X[i - 1] + h / 2, Y[i - 1] + k1 / 2);
		k3 = h * function(X[i - 1] + h / 2, Y[i - 1] + k2 / 2);
		k4 = h * function(X[i - 1] + h, Y[i - 1] + k3);
		dY[i - 1] = k1 / 6 + k2 / 3 + k3 / 3 + k4 / 6;
		Y[i] = Y[i - 1] + dY[i - 1];
		printf_s("i = %d, X = %lf, Y = %lf\n", i, X[i], Y[i]);
	}
	//printf_s("Answer: %lf\n", function(X[10], Y[10]));
}

void AdamsMeth(double a, double b, double y0, double h, double exp) {
	double k[4], y[4], x[4];
	double X = a, Y = y0, dy, new_y;
	int count = 1;
	//Palkas();
	printf_s("\nAdams method:\n%d. x = %.2lf\ty = %.2lf\n", 0, X, Y);
	for (int i = 0; i < 4; i++) {
		y[i] = Y;
		x[i] = X;
		k[0] = h * function(X, Y);
		k[1] = h * function(X + h / 2, Y + k[0] / 2);
		k[2] = h * function(X + h / 2, Y + k[1] / 2);
		k[3] = h * function(X + h, Y + k[2]);
		dy = (k[0] + 2 * k[1] + 2 * k[2] + k[3]) / 6;
		Y += dy;
		X += h;
		printf("%d. x = %.2lf\ty = %.2lf\n", count, X, Y);
		count++;
	}
	while (X <= b) {
		Y = y[3] + h / 24 * (function(x[3], y[3]) * 55 - function(x[2], y[2]) * 59 + function(x[1], y[1]) * 37 - function(x[0], y[0]) * 9);
		X += h;
		while (1) {
			new_y = y[3] + h / 24 * (9 * function(X, Y) + 19 * function(x[3], y[3]) - 5 * function(x[2], y[2]) + function(x[1], y[1]));
			if (fabs(function(X, new_y) - function(X, Y)) < exp) break;
			Y = new_y;
		}
		Y = y[3] + h / 24 * (9 * function(X, new_y) + 19 * function(x[3], y[3]) - 5 * function(x[2], y[2]) + function(x[1], y[1]));
		printf("%d. x = %.2lf\ty = %.2lf\n", count, X, Y);
		count++;
		for (int i = 0; i < 3; i++) {
			y[i] = y[i + 1];
			x[i] = x[i + 1];
		}
		y[3] = Y;
		x[3] = X;
	}
}