function [w,t] = AdamsAdaptPCv(FunFcnIn, Intv, alpha, tol,stepsize)
% On input: 
%   FunFcnIn is the name of function to be integrated
%   interv is the interval to be integrated over
% 
%   The problem: y' = f(t,y), y(a) = alpha, a<= t <= b
%   where Intv = [a b], and a call to FunFcnIn with 
%   argument (t, y) returns f(t,y).
%
%   AdamsAdaptPC uses the Adams Variable Step-size PC method to solve
%   the above problem, to a given tolerance. 
%   hmin and hmax are the minimum and maximum step sizes.
%
% This is the vector version of AdamsAdaptPC.
%
% On output
%   t contains the (unequi-spaced) mesh points, w the
%   function values at these points. 
%
% Written by Ming Gu for Math 128A, Fall 2008
% 
[FunFcn,msg] = fcnchk(FunFcnIn,0);
if ~isempty(msg)
    error('InvalidFUN',msg);
end
flg  = 0;
a    = Intv(1);
b    = Intv(2);
hmin = stepsize(1);
hmax = stepsize(2);
h    = min(hmax,(b-a)/4);
if (h<=0)
    msg = ['Illegal hmax or interval'];
    error('InvalidFUN',msg);
end
p = [-9/24  37/24 -59/24 55/24]';
c = [ 1/24 -5/24   19/24 9/24 ]';
%
% RK4v for the first 3 steps
%
nflg      = 1;
[f, w, t] = RKv(FunFcnIn, a,h,alpha);
%
% main loop
%
flg = 0;
i   = 4;
last= 0;
while (flg == 0)
    wi   = w(:,i);
    ti   = t(i);
    tnew = ti + h;
    wp   = wi + h*(f(:,i-3:i)*p);
    fp   = FunFcn(tnew,wp);
    wc   = wi + h*([f(:,i-2:i),fp]*c);
    sigma= max(eps,19*norm(wc-wp)/(270*h));
    if (sigma <= tol)
        i    = i + 1;
        t(i) = tnew;
        w(:,i) = wc;
        f(:,i) = FunFcn(tnew,wc);
        if (last == 1)
            return;
        end
        nflg = 0;
        if (sigma <= 0.1*tol | tnew+h>b)
            q  = (tol/(2*sigma))^(1/4);
            h  = min(hmax,min(q,4)*h);
            if (tnew + 4 * h > b)
                h = (b-tnew)/4;
                last = 1;
            end
            nflg = 1;
            [f1, w1, t1] = RKv(FunFcnIn, t(i), h, w(:,i));
            f = [f(:,1:i),f1(:,2:end)];
            w = [w(:,1:i),w1(:,2:end)];
            t = [t(1:i);t1(2:end)];
            i = i + 3;
        end
    else
        q  = (tol/(2*sigma))^(1/4);
        h  = max(q,0.1)*h;
        if (h < hmin)
            disp(['hmin exceeded']);
            flg = 1;
            return;
        end
        if (nflg == 1)
            i = i - 3;
        end
        nflg = 1;
        [f1, w1, t1] = RKv(FunFcnIn, t(i), h, w(:,i));
        f = [f(:,1:i),f1(:,2:end)];
        w = [w(:,1:i),w1(:,2:end)];
        t = [t(1:i);t1(2:end)];
        i = i + 3;
    end
end

function [f,w, t] = RKv(FunFcnIn, a,h,alpha)
%
% performs 3 steps of RK4v
%
[FunFcn,msg] = fcnchk(FunFcnIn,0);
if ~isempty(msg)
    error('InvalidFUN',msg);
end
m    = length(alpha);
w    = zeros(m,4);
f    = zeros(m,4);
w(:,1) = alpha;
h2   = h/2;
t    = a+(0:3)'*h;
for i = 1:3
    f(:,i) = FunFcn(t(i),w(:,i));
    k1 = h* f(:,i);
    k2 = h* FunFcn(t(i)+h2,w(:,i)+k1/2);
    k3 = h* FunFcn(t(i)+h2,w(:,i)+k2/2);
    k4 = h* FunFcn(t(i)+h,w(:,i)+k3);
    w(:,i+1) = w(:,i) + (k1+2*k2+2*k3+k4)/6;
end
f(:,4) = FunFcn(t(4),w(:,4));