Matlab example 2: Explicit vs Implicit 2 steps Adams methods

clear all
close all

f = @(x,y) y^2;
real_y = @(x) 1/(1-x);

I = [0 0.99];
numsteps = 20;
h = (I(2)-I(1))/numsteps;

x(1) = 0;
y_e(1) = 1;
y_i(1) = 1;
y(1) = 1;

figure
plot(x,y_e,'-*b');
hold on
plot(x,y_i,'-*r');
plot(x,y,'--k');
xlim(I)
ylim([0, 20])
axis square

F_e = [f(x(1),y_e(1))];
F_i = [f(x(1),y_i(1))];

for i = 1 : numsteps-1
    % Compute the new grid point
    x(i+1) = x(i) + h;
    % Compute the new approximation with EAM
    y_e(i+1) = EAM2(y_e(i),F_e,h);
    F_e = [f(x(i+1),y_e(i+1)) F_e];
    % Compute the new approximation with IAM
    y_i(i+1) = PEC_IAM2(x(i+1),y_i(i),F_i,h,f);
    F_i = [f(x(i+1),y_i(i+1)) F_i];
    % Compute the value of the true solution
    y(i+1) = real_y(x(i+1));

    plot(x,y_e,'-*b');
    plot(x,y_i,'-*r');
    plot(x,y,'--k');
    xlim(I)
    ylim([0, 20])
    axis square
    pause(.1)
end
legend('EAM(2)', 'IAM(2) with PEC', 'True Solution', 'location','northwest');

hold off






function y = EAM2(yn, F, h)

    % INPUT:   yn : solution approx on x_n
    %          F  : vector of previous values F = [f_n, f_{n-1}]
    %          h  : step size
    % OUTPUT:  y  : new value y_{n+1}

    if length(F) >= 2
        y = yn + h * .5 * (3*F(1) - F(2));
    elseif length(F) == 1
        y = yn + h * F(1); %Explicit Euler method
    end    

end

function y = PEC_IAM2(xn1, yn, F, h, f)

    % INPUT:   xn1 : next grid point x_{n+1}
    %          yn  : solution approx on x_n
    %          F   : vector of previous values F = [f_n, f_{n-1}]
    %          h   : step size
    %          f   : function handler for f(x,y)
    % OUTPUT:  y   : new value y_{n+1}

    % Predictor
    yn1 = EAM2(yn, F, h);

    % Evaluation
    F = [f(xn1,yn1) F];

    % Corrector (IAM2)
    y = yn + h * .5 * (F(1) + F(2));

end