function[x,u] = AdamsMoulton(a,b,u0,N)

h = (b-a)/N;

[xrk,urk] = rk(a,a+3*h,u0,3);
for i = 1:4
    x(i) = xrk(i);
    u(i) = urk(i);
end

for i = 5:N+1
    x(i) = a+(i-1)*h;
    u(i) = u(i-1) + h/24*( 55*fw(x(i-1),u(i-1)) ...
                         - 59*fw(x(i-2),u(i-2)) ...
                         + 37*fw(x(i-3),u(i-3)) ...
                         -  9*fw(x(i-4),u(i-4)) );
    u(i) = u(i-1) + h/24*(  9*fw(x(i),u(i)) ...
                         + 19*fw(x(i-1),u(i-1)) ...
                         -  5*fw(x(i-2),u(i-2)) ...
                         +    fw(x(i-3),u(i-3)) );                
end