As follows
f=@(t,y) 3*t+y/t;
alpha=5;
a=1;
b=2;
n=3;
[t, w, h] = abs2(f, a, b, alpha, n);
plot(t,w,'-o'),grid
xlabel('t'),ylabel('w')
disp(['h = ' num2str(h)])
function [t, w, h] = abs2(f, a, b, alpha, n)
%AB2 Two-step Adams Bashforth method
% [t, w, h] = ab2(f, a, b, alpha, n) performs the two-step Adams Bashforth
% method for solving the IVP y' = f(t,y) with initial condition y(a) = alpha
% taking n steps from t = a to t = b. The first step from t = a to t = a + h
% is performed using the modified Euler method.
h=(b-a)/n;
t=a:h:b;
w=zeros(1,length(t));
w(1)=alpha;
%modified euler method
k1=h*f(t(1),w(1));
k2=h*f(t(1)+h,w(1)+k1);
w(2)=w(1)+1/2*(k1+k2);
for i=2:length(t)-1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
w(i+1)=w(i)+(3/2)*h*f(t(i),w(i))-.5*h*f(t(i-1),w(i-1));
end
end
