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