Variable step size Implicit Adams method

Using the recursions for the coefficients $\Phi_j(n)$, $\Psi_j(n)$ and $g_j(n)$, which we recall below

$\Phi_0(n)=\Psi_0(n)=f_n$

$\Phi_{j+1}(n) = \Phi_j(n)-\Psi_j(n-1)$

$\Psi_{j+1}(n) = \beta_{j+1}(n)\Phi_{j+1}(n)$

$\beta_{j}(n) = \beta_{j-1}(n)\Big(\frac{x_{n+1}-x_{n-j+1}}{x_n-x_{n-j}}\Big),\qquad$ $\beta_0(n)=1$

$Q(0,q)=1/q, \qquad$ $Q(1,q)=\frac{1}{q(q+1)}$

$Q(j,q) = Q(j-1,q)-Q(j-1,q+1)$

$g_j(n)= Q(j,1)$

Implement a 1-step and 2-step IAM with adaptive step size to solve the ODE $y'=y^2$, $y(0)=1$, using the scheme

$y_{n+1} = p_{n+1} + h_n g_k(n) \Phi_k(n+1),\qquad$ $p_{n+1} = y_n + h_n \sum_{j=0}^{k-1}g_j(n)\Psi_j(n)$

Matlab DEMO

clear all

%  F = history of function evaluations
%  X = history of points
%  Y = history of solution approximations

f = @(x,y) y^2;
sol = @(x) 1./(1-x); %true solution
X = [0.01 0.02 0.03];
Y = [sol(X(1)) sol(X(2)) sol(X(3))];
F = [f(X(1),Y(1)) f(X(2),Y(2)) f(X(3),Y(3))];

hnew = 0.13;
h = hnew;
tol = 1e-5;
%  Interval [0,1]

while X(end) < 0.99

    [y,hnew] = corrector(hnew,F,X,Y,f,tol);

    X = [X X(end) + hnew];
    Y = [Y y];
    F = [F f(X(end),Y(end))];
    h = [h hnew];
end

Y = Y(1:end-1);
X = X(1:end-1);

close all
subplot(121)
plot(X,Y,'-or');
hold on
plot(X,sol(X),'--k', 'linewidth',3)
title('Approx solution')
subplot(122)
semilogx(h,'-')
title('Sequence of steps h_n')


function [y,hnew] = corrector(hnew,F,X,Y,f,tol)

    for it = 1:100
        [p,PSI,G] = predictor(hnew,F,X,Y);
        Xn1 = X(end)+hnew;
        G(3) = 1/2 - (1/6)*(hnew/(Xn1-X(end-1)));
        PHI3 = f(Xn1,p) - F(end) - PSI(2);
        y = p + hnew * G(3) * PHI3;

        hn = hnew;

        %%% We compute an approximation LE for the local error at the next point
        %%% and compute a new step length as a consequence
        % compute PHI_3(n+1)
        b1n = hn / (X(end)-X(end-1));
        c1 = (Xn1 - X(end-1))/(X(end)-X(end-2));
        c2 = (X(end)-X(end-1))/(X(end-1)-X(end-2));
        PHI4 = PHI3 - b1n * c1 * (F(end)-F(end-1)-c2*(F(end-1)-F(end-2)));
        % compute g_3(n)
        G(4) = G(3) - (hn/(Xn1-X(end-2))) * (1/6 - 1/12 * (hn / (Xn1 - X(end-1))));
        % compute LE
        LE = hn * abs(PHI4 * (G(4)-G(3)) );
        % define a new hnew
        hnew = hn * (tol/LE)^(1/4);

        %%% If the approximate LE is smaller than the tollerance we accept
        if LE < tol
            break;
        end

    end
end

function [p,PSI,G] = predictor(hnew,F,X,Y)
    G(1) = 1;
    G(2) = 1/2;
    PSI(1) = F(end);
    Xn1 = X(end) + hnew;
    PSI(2) = (F(end)-F(end-1))*((Xn1-X(end))/(X(end)-X(end-1)));
    p = Y(end) + hnew * (G(1)*PSI(1) + G(2)*PSI(2));
end