function [t,w] = mvsspc(f,a,alpha,h0,tol,N)
%Implements a simplified variable stepsize predictor corrector method using four step
%Adams-Bashform (explicit) and three step Adams-Moulton (implicit). 




t = zeros(1,N);
w = zeros(length(alpha),N);
t(1) = a;
w(:,1) = alpha;
h = h0;
changestep = 1; %1 indicates we changed the stepsize, 0 if not.
i = 1; %i indicates that we've computed t and w up to component i

while (i < N)
    
    if changestep == 1
        [t(i:i+3), w(:,i:i+3)] = myrk4(f,t(i),t(i)+3*h,w(:,i),3);
        i = i+3;
    end
    
    WP = w(:,i) + h/24*(55*f(t(i),w(:,i)) - 59*f(t(i-1),w(:,i-1)) + ...
        37*f(t(i-2),w(:,i-2)) - 9*f(t(i-3),w(:,i-3)) );
    WC = w(:,i) + h/24*(9*f(t(i)+h,WP) + 19*f(t(i),w(:,i)) - ...
        5*f(t(i-1),w(:,i-1)) + f(t(i-2),w(:,i-2)) );
    sigma = 19*norm( WC - WP,inf )/(270*h);
    
    if sigma <= tol %Result accepted
        w(:,i+1) = WC;
        t(i+1) = t(i) + h;
        i = i + 1;
        changestep = 0;
        if sigma <= .1*tol %Result too accurate, increase h.
            h = h*(tol/(2*sigma))^(1/4);
            changestep = 1;
        end
    else %Result rejected
        if changestep == 1 
            i  = i - 3;
        end
        h = h*(tol/(2*sigma))^(1/4);
        changestep = 1;
    end
end
end