% dy/dx = x +y
% dy/dx = x.^2 -y

% Zawsze musimy miec rownanie w postaci normalnej (dy/dx=... / y'=...)

syms y(x)
fA = dsolve(diff(y) == (x-y)/2, y(0) == 1)
fn1 = @(x, y) (x + y)./2
d = @(x) 3* exp(-x./2) - 2

clf;

a = 0
b = 3
f = fn1
tspan = [a b]
eps = 10^7

n = (b - a)/0.125
[x1, y1] = abm3(f, tspan, 1, n, eps)
% [x2, y2] = ab3(f, tspan, 1, eps)
[x3, y3] = ode45(f, tspan, 1)

title("Krok 0.125")
subplot(2,1,1); 

hold on
plot(x1, y1)
% plot(x2, y2)
plot(x3, y3)
fplot(fA)
legend('abm3', 'ode45', 'dokl')

subplot(2,1,2);
hold on

plot(x1, abs(y1 - f(x1)))
plot(x3, abs(y1 - f(x3)))

hold off

function [x, y] = abm3(f, tspan, y0, n, eps)
    a = tspan(1);
    b = tspan(2);
    h = (b - a)/n;
    x = (a : h : b);

    % Euler
    y(1) = y0;
    for i = 1 : 2
        k1 = h * feval(f, x(i), y(i));
        k2 = h * feval(f, x(i) + h/2, y(i) + k1/2);
        y(i+1) = y(i) + k2;
    end
    z0 = feval(f, a, y0);
    z(1) = feval(f, x(1), y(1));
    z(2) = feval(f, x(2), y(2));

    % Predyktor      
    y(3) = y(2) + h / 12 * (23 * z(2) - 16 * z(1) + 5 * z0);
    for i = 3 : n
        z(i) = feval(f, x(i), y(i));
    
        yy = [y(i) + h / 12 * (23 * z(i) - 16 * z(i-1) + 5 * z(i-2))];
        s = 1;

        for awaryjna = 1:100
            a = feval(f, x(i), yy(s));
            b = feval(f, x(i - 1), y(i - 1));
            c = feval(f, x(i - 2), y(i - 2));
            % Korektor
            yy_c = y(i-1) + 1/12 * h * (5 * a + 8 * b - c);
            yy = [yy yy_c];

            if abs(yy(length(yy)) - yy(length(yy - 1))) < eps
                break;
            end
            s = s + 1
        end
        y(i+1) = yy(length(yy));
    end