Create residual viscosity package +rv and generalize the ResidualViscosity class - Generalize residual viscosity, by passing user-defined flux and calculating the time derivative outside of the update. - Create separate RungekuttaRV specifically using interior RV updates - Separate the artifical dissipation operator from the scheme AdvectionRV1D so that the same scheme can be reused for creating the diff op used by the ResidualViscosity class

% Subs for a symfun
% f remains a symbolic function. If any of it's arguments is eliminated
% it is removed from the argument list while preserving the order of the
% other arguments
function f = subsSymfun(f, old, new)
    args = argnames(f);

    newExpr = subs(f, old, new);
    vars = symvar(subs(args, old, new));

    newArgs = args(ismember(args,vars));

    f = symfun(newExpr, newArgs);
end