function [t,w]=ABM(y0,a,b,h)
% ABM(1,0,1,0.1);

N=(b-a)/h;
t=linspace(a,b,N+1);
w(1)=y0;
for i=1:3
    k1=h*defunc(t(i),w(i));
    k2=h*defunc(t(i)+h/2,w(i)+k1/2);
    k3=h*defunc(t(i)+h/2,w(i)+k2/2);
    k4=h*defunc(t(i)+h,w(i)+k3);
    w(i+1)=w(i)+(k1+2*k2+2*k3+k4)/6;
end
for i=4:N
   w(i+1)=w(i)+h/24*(55*defunc(t(i),w(i))-59*defunc(t(i-1),w(i-1))+37*defunc(t(i-2),w(i-2))-9*defunc(t(i-3),w(i-3)));
   w(i+1)=w(i)+h/24*(9*defunc(t(i+1),w(i+1))+19*defunc(t(i),w(i))-5*defunc(t(i-1),w(i-1))+defunc(t(i-2),w(i-2)));
end
[T,Y] = ODE45('defunc',t,y0);
plot(T,Y)
hold on
plot(t,w,'o')
hold off
error=max(abs(w-Y'));
