AdamsMoulton:=proc(k)
local p,i,dt,t;
  p:=interp([seq(i*dt,i=0..k)],[seq(f[i],i=0..k)],t);
  y[k]-y[k-1]=factor(int(p,t=(k-1)*dt..k*dt));
end: