restart;
# Find Adams-Bashforth method of order s
s:=4;
t:=n->t_0+n*h;
L:=k->product((tau-t(n-j))/(t(n-k)-t(n-j)),j=0..k-1)
     *product((tau-t(n-j))/(t(n-k)-t(n-j)),j=k+1..s-1);
p:=0:
for i from 0 to s-1
do
    p:=p+simplify(L(i))*f[n-i];
od:
p:=collect(p,tau);
AB:=simplify(int(p,tau=t(n)..t(n+1)));
y[n+1]=y[n]+AB;