function oscillator1d_avi_3pot(ts1,ts2,ts3,k1,k2,k3,x0,v0,gamma,T)

% This program computes AVI propagator matrix for the system 
%     a + gamma v + (k1+k2+k3) x = 0
% and advances the position and velocity from t=0 to t=T.

% inputs:
% ts1 - small timestep
% ts2 - medium timestep
% ts3 - large timestep
% k1 - stiff spring constant
% k2 - medium spring constant
% k3 - soft spring constant
% T - total simulation time
% x0 - initial position
% v0 - initial velocity
% gamma - friction coefficient

% error checking
if (ts1 <= 0 || ts2 <= 0 || ts3 <= 0 || ts1 > ts2 || ts2 > ts3)
    disp('Must have 0 < ts1 <= ts2 <= ts3');
    return;
end
if (k1 <= 0 || k2 <= 0 || k3 <= 0 || k1 < k2 || k2 < k3)
    disp('Must have k1 >= k2 >= k3 > 0');
    return;
end
if (T <= 0)
    disp('Must have total simulation time T > 0');   
    return;
end

% initialization
i = 1;              % index for t and xv
t(i) = 0;           % time vector
xv(:,1) = [x0;v0];  % phase space vector
tslow = 0;          % time of last update of slow step
tmed  = 0;          % time of lase update of med step
tfast = 0;          % time of last update of fast step

% sets up priority queue
pq(1,:) = [ts1 1];
pq(2,:) = [ts2 2];
pq(3,:) = [ts3 3];

% sets up the matrices that are time-invariant
w = (1 - ts1/2*gamma)/(1 + ts1/2*gamma);      % extra factor for friction
Vfast = [1 0; -ts1/2*k1*(1+w) w];         % velocity update matrices
Vmed  = [1 0; -ts2*k2 1];
Vslow = [1 0; -ts3*k3 1];

% advance velocity by half time step for each potential
xv(:,1) = [1 0; -0.5*ts1*k1/2*(1+w) w]*xv(:,1);
xv(:,1) = [1 0; -0.5*ts2*k2 1]*xv(:,1);
xv(:,1) = [1 0; -0.5*ts3*k3 1]*xv(:,1);

% integrate until time T (1.e-10 is there to account for roundoff errors)
while (tfast < T - 1.e-10 || tmed < T -1.e-10 || tslow < T - 1.e-10)
    if (pq(1,2) == 1)                    % update fast potential
        dt = tfast+ts1-t(i);             % time interval
        Xfast = [1 dt; 0 1];         
        xv(:,i+1) = Vfast*Xfast*xv(:,i); % update phase space vector
        tfast = tfast + ts1;             % update time of potential
        t(i+1) = tfast;                  % advances current time of particle
        
        pq(1,1) = pq(1,1) + ts1;         % schedule next update
    elseif (pq(1,2) == 2)                % update med potential
        dt = tmed+ts2-t(i);              % time interval
        Xmed = [1 dt; 0 1];         
        xv(:,i+1) = Vmed*Xmed*xv(:,i);   % update phase space vector
        tmed = tmed + ts2;               % update time of potential
        t(i+1) = tmed;                   % advances current time of particle
        
        pq(1,1) = pq(1,1) + ts2;         % schedule next update
    else                                 % update slow potential
        dt = tslow+ts3-t(i);             % time interval
        Xslow = [1 dt; 0 1];         
        xv(:,i+1) = Vslow*Xslow*xv(:,i); % update phase space vector
        tslow = tslow + ts3;             % update time of potential
        t(i+1) = tslow;                  % advances current time of particle
        
        pq(1,1) = pq(1,1) + ts3;         % schedule next update
    end
    
    % increment the index
    i = i+1;
    
    % sort priority queue as needed
    if (pq(1,1) > pq(3,1))
        tmp = pq(1,:);
        pq(1:2,:) = pq(2:3,:);
        pq(3,:) = tmp;
    elseif (pq(1,1) > pq(2,1))
        tmp = pq(1,:);
        pq(1,:) = pq(2,:);
        pq(2,:) = tmp;
    end
end

% plots the position of the oscillator
figure(1); clf
plot(t,xv(1,:));
xlabel('Time');
ylabel('Position');