function march( c2R, c2I, rho, x, Ik, deltak, rr )

% Solve system for a sequence of k-values

global NMedia N H
global Pulse fMin fMax tStart alpha beta V
global tout deltat
global cLow cHigh
global omega Bdry
global WS Isd WR Ird Z
global Green RTSrd RTSrr

MaxIT   = 1e6;
deltat2 = deltat^2;

% March a single spectral component forward in time
NTot1 = sum( N ) + 1;
AD0   = zeros( NTot1, 1 );
AE0   = zeros( NTot1, 1 );
AD1   = zeros( NTot1, 1 );
AE1   = zeros( NTot1, 1 );
AD2   = zeros( NTot1, 1 );
AE2   = zeros( NTot1, 1 );
Ke    = zeros( 2, 2 );
Ce    = zeros( 2, 2 );
Me    = zeros( 2, 2 );

% Assemble mass and stiffness matrices

rkT  = sqrt( x );

Node = 1;
L    = 1;

for medium = 1 : NMedia
    hElt = H( medium );

    for I = 1 : N( medium )
        rhoElt = ( rho( L ) + rho( L + 1 ) ) / 2.0;
        c2RElt = ( c2R( L ) + c2R( L + 1 ) ) / 2.0;
        c2IElt = ( c2I( L ) + c2I( L + 1 ) ) / 2.0;
        rhoh   = rhoElt * hElt;

        % Elemental stiffness matrix ( +k2 term )
        y = 1.0 + i * rkT * V * c2IElt;
        z = hElt * x * ( y - V * V / c2RElt ) / rhoElt / 6.0;
        Ke( 1, 1 ) =  y / rhoh + (3.0-alpha) * z;
        Ke( 1, 2 ) = -y / rhoh +      alpha  * z;
        Ke( 2, 2 ) =  y / rhoh + (3.0-alpha) * z;

        % Elemental damping matrix
        xx = hElt * x / rhoElt / 6.0;
        z  = hElt * 2.0 * i * rkT * V / c2RElt / rhoElt / 6.0;
        Ce( 1, 1 ) = c2IElt * ( 1.0/rhoh + (3.0-alpha) * xx) + (3.0-alpha) * z;
        Ce( 1, 2 ) = c2IElt * (-1.0/rhoh +      alpha  * xx) +      alpha  * z;
        Ce( 2, 2 ) = c2IElt * ( 1.0/rhoh + (3.0-alpha) * xx) + (3.0-alpha) * z;

        % Elemental mass matrix
        Me( 1, 1 ) = (3.0-alpha) * hElt / ( rhoElt * c2RElt ) / 6.0;
        Me( 1, 2 ) =      alpha  * hElt / ( rhoElt * c2RElt ) / 6.0;
        Me( 2, 2 ) = (3.0-alpha) * hElt / ( rhoElt * c2RElt ) / 6.0;

        % A2 matrix
        AD2( Node   ) = AD2( Node ) + ...
            Me(1,1) + 0.5 * deltat * Ce(1,1) + beta * deltat2 * Ke(1,1);
        AE2( Node+1 ) = Me(1,2) + 0.5 * deltat * Ce(1,2) + beta * deltat2 * Ke(1,2);
        AD2( Node+1 ) = Me(2,2) + 0.5 * deltat * Ce(2,2) + beta * deltat2 * Ke(2,2);

        % A1 matrix
        AD1( Node   ) = AD1( Node ) + ...
            2.0*Me(1,1) - ( 1.0 - 2.0*beta ) * deltat2 * Ke(1,1);
        AE1( Node+1 ) = 2.0*Me(1,2) - ( 1.0 - 2.0*beta ) * deltat2 * Ke(1,2);
        AD1( Node+1 ) = 2.0*Me(2,2) - ( 1.0 - 2.0*beta ) * deltat2 * Ke(2,2);

        % A0 matrix
        AD0( Node   ) = AD0( Node ) ...
            -Me(1,1) + 0.5 * deltat * Ce(1,1) - beta * deltat2 * Ke(1,1);
        AE0( Node+1 ) = -Me(1,2) + 0.5 * deltat * Ce(1,2) - beta * deltat2 * Ke(1,2);
        AD0( Node+1 ) = -Me(2,2) + 0.5 * deltat * Ce(2,2) - beta * deltat2 * Ke(2,2);

        Node = Node + 1;
        L    = L    + 1;
    end
    L = L + 1;
end

%factor( NTot1, AD2, AE2, RV1, RV2, RV3, RV4 )  % * Factor A2 *

% * Initialize pressure vectors *
U0 = zeros( NTot1, 1 );
U1 = zeros( NTot1, 1 );

% Initializate parameters for bandpass filtering of the source time series
if ( Pulse(4:4) == 'L' || Pulse(4:4) == 'B' )
    fLoCut = rkT * cLow / ( 2.0 * pi );
else
    fLoCut = 0.0;
end

if ( Pulse(4:4) == 'H' || Pulse(4:4) == 'B' )
    fHiCut = rkT * cHigh / ( 2.0 * pi );
else
    fHiCut = 10.0 * fMax;
end

time   = 0.0;
IniFlag = 1;   % need to tell source routine to refilter for each new value of k

% Forward, march
n  = length( AE0 );
AF0 = [ AE0( 2 : end ); 0 ];    % to conform to Matlab convention relating diagonals to the matrix
AF1 = [ AE1( 2 : end ); 0 ];
A0 = spdiags( [ AF0 AD0 AE0 ], -1:1, n, n );
A1 = spdiags( [ AF1 AD1 AE1 ], -1:1, n, n );
Itout = 1;

for Itime = 1 : MaxIT
    time  = tStart + ( Itime - 1 ) * deltat;
    NTot1 = length( U0 );

    % ** Form U2TEMP = A1*U1 + A0*U0 + S **
    U2 = A1 * U1 + A0 * U0;

    %Ntpts = 1;
    %timeV( 1 ) = time; 	% source expects a vector, not a scalar
    %source( timeV, ST, sd, Nsd, Ntpts, omega, fLoCut, fHiCut, Pulse, PulseTitle, IniFlag );
    [ ST, PulseTitle ] = cans( time, omega, Pulse );
    
    ST = deltat2 * ST * exp( -i * rkT * V * time );
    js = Isd;    % assumes source in first medium
    U2( js   ) = U2( js   ) + ( 1.0 - WS ) * ST;
    U2( js+1 ) = U2( js+1 ) +         WS   * ST;

    % Solve A2*U2 = U2TEMP
    if ( alpha == 0.0 && beta == 0.0 )   % explicit solver
        U2 = U2 ./ AD2;
    else
        backsb( NTot1, RV1, RV2, RV3, RV4, U2 );   % implicit solver
    end
    
    % Do a roll
    U0 = U1;
    U1 = U2;

    % Boundary conditions (natural is natural)
    if ( Bdry.Top.Opt(2:2) == 'V' )
        U1( 1 )     = 0.0;
    end
    if ( Bdry.Bot.Opt(1:1) == 'V' )
        U1( NTot1 ) = 0.0;
    end

    % Extract solution (snapshot, vertical or horizontal array)
    while ( ( Itout <= length( tout ) ) && ( time + deltat >= tout( Itout ) ) )
        % above can access tout( ntout + 1 ) which is outside array dimension but result is irrelevant
        wt = ( tout( Itout ) - time ) / deltat;  % Weight for temporal interpolation:

        switch ( Bdry.Top.Opt(4:4) )
            case ( 'S' ) % snapshot
                const = exp( i * rkT * V * tout( Itout ) );
                UT1 = U0( Ird ) + WR .* ( U0( Ird+1 ) - U0( Ird ) );   % Linear interpolation in depth
                UT2 = U1( Ird ) + WR .* ( U1( Ird+1 ) - U1( Ird ) );
                Green( Itout, :, Ik ) = const * ( UT1 + wt * ( UT2 - UT1 ) );   % Linear interpolation in time

            case ( 'D' ) % RTS (vertical array)
                const =  sqrt( 2.0 ) * deltak * sqrt( rkT ) * exp( i*( rkT * ( V * tout( Itout ) - rr( 1 ) ) + pi / 4.0 ) );
                UT1 = U0( Ird ) + WR .* ( U0( Ird+1 ) - U0( Ird ) );   % Linear interpolation in depth
                UT2 = U1( Ird ) + WR .* ( U1( Ird+1 ) - U1( Ird ) );
                U   = UT1 + wt * ( UT2 - UT1 );                        % Linear interpolation in time
                RTSrd( :, Itout ) = RTSrd( :, Itout ) + real( const * U );

            case ( 'R' ) % RTS (horizontal array)
                if ( time >= tout( 1 ) )
                    % Linear interpolation in depth
                    I     = Ird( 1 );
                    T     = ( rd( 1 ) - Z( I ) ) / ( Z( I + 1 ) - Z( I ) );
                    UT1   = U0( I ) + T * ( U0( I+1 ) - U0( I ) );
                    UT2   = U1( I ) + T * ( U1( I+1 ) - U1( I ) );
                    U     = UT1 + wt * ( UT2 - UT1 );                       % Linear interpolation in time
                    const = sqrt( 2.0 ) * exp( i * rkT * V * time );
                    RTSrr( :, Itout ) = RTSrr( :, Itout ) + real( deltak * U *  ...
                        const * exp( i * ( -rkT * r + pi / 4.0 ) ) * sqrt( rkT / rr ) );
                end
        end
        Itout = Itout + 1;
    end
    if ( Itout > length( tout ) )    % jump out of while loop and do next wavenumber
        break;
    end
end