% Takes one time step of size dt using the rungekutta method
% starting from v and where the function F(v,t,RV) gives the
% time derivatives. coeffs is a struct holding the RK coefficients
% for the specific method. RV is the residual viscosity which is updated
% in between the stages and after the updated solution is computed.
function v = rungekuttaRV(v, t , dt, F, RV, DvDt, coeffs)
    % Move one stage outside to avoid branching for updating the
    % residual inside the loop.
    k = zeros(length(v), coeffs.s);
    k(:,1) = F(v,t,RV.evaluate(v,DvDt(v)));

    % Compute the intermediate stages k
    for i = 2:coeffs.s
        u = v;
        for j = 1:i-1
            u = u + dt*coeffs.a(i,j)*k(:,j);
        end
        %RV.update(0.5*(u+v),(u-v)/(coeffs.c(i)*dt)); % Crank-Nicholson for time discretization
        k(:,i) = F(u,t+coeffs.c(i)*dt, RV.evaluate(u,DvDt(u)));
    end

    % Compute the updated solution as a linear combination
    % of the intermediate stages.
    u = v;
    for i = 1:coeffs.s
        u = u + dt*coeffs.b(i)*k(:,i);
    end
    %RV.update(0.5*(u+v),(u-v)/dt); % Crank-Nicholson for time discretization
    v = u;
end