function [x,y] = ab3(f, tspan, y0, n)
    a = tspan(1);
    b = tspan(2);

    h = (b - a) / n;
    x = (a : h : b);

    y(1) = y0;
    for i = 1 : 2
        k_1 = h * feval(f, x(i), y(i));
        k_2 = h * feval(f, x(i) + h / 2, y(i) + k_1/2);

        y(i+1) = y(i) + k_2;
    end

    z0 = feval(f, a, y0);
    z(1) = feval(f, x(1), y(1));
    z(2) = feval(f, x(2), y(2));

    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));
        y(i+1) = y(i) + h / 12 * (23 * z(i) - 16 * z(i-1) + 5 * z(i-2));
    end

function [t, y] = euler(f, tspan, y0, N)
    yy = [y0];
    xx = [0];

    a = tspan(1);
    b = tspan(2);
    h = (b - a) / N;

    for i = 2:N
        current_x = xx(i-1) + h;
        current_y = yy(i-1) + h * f(current_x, yy(i-1));
        xx = [xx current_x];
        yy = [yy current_y];
        if current_x > 20
            break;
        end;
    end
    t = xx;
    y = yy;

function [t, y] = euler_m(f, tspan, y0, N)
    yy = [y0];
    xx = [0];

    a = tspan(1);
    b = tspan(2);
    h = (b-a)/N;
    for i = 2:N
        current_x = xx(i-1) + h;
        current_y = yy(i-1) + h * f(current_x, yy(i-1));
        xx = [xx current_x];
        yy = [yy current_y];
        if current_x > 20
            break;
        end;
    end
    t = xx;
    y = yy;

syms y(x)
dsolve(diff(y) == x+y, y(0) == 1)

dsolve(diff(y) == x^2 - y, y(0) == 1)

fn1 = @(x, y) x+y
fn2 = @(x, y) x^2 * x - y * x
fd1 = @(x) 2 * exp(x) - x -1 
fd2 = @(x) x^2 -exp(-x) - 2 * x +x
a = 0
b = 10
tspan = [0 10];

[x1a, y1a] = euler(fn1, tspan, 1, (b - a) / 0.01)
[x1b, y1b] = euler(fn2, tspan, 1, (b - a) / 0.01)

[x2a, y2a] = euler_m(fn1, tspan, 1, (b - a) / 0.01)
[x2b, y2b] = euler_m(fn2, tspan, 1, (b - a) / 0.01)

[x3a, y3a] = ab3(fn1, tspan, 1, (b-a) / 0.01)
[x3b, y3b] = ab3(fn2, tspan, 1, (b-a) / 0.01)

hold on
plot(x1a, y1a)
plot(x2a, y2a)
plot(x3a, y3a)
plot(x4a, y4a)

hold off

hold on
plot(x1b, y1b)
plot(x2b, y2b)
plot(x3b, y3b)
plot(x4b, y4b)

hold off

plot(x1a, y1a)

err1 = abs(y4a - y1a)
err2 = abs(y4a - y1a)
err3 = abs(y4a - y1a)