function ts = getRvTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0)
    switch opt.method
        case 'rkRv'
            ts = rkRvTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0);
        case 'rkRvBdf'
            ts = rkRvBdfTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0);
        case 'rkRvMs'
            ts = rkRvMsTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0);
        case 'rkRvMg'
            ts = rkRvMgTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0);
        case 'rkRvInstage'
            ts = rkRvInstageTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0);
        otherwise
            error('Timestepping method ''%s'' not supported',method);
    end
end

function fhData = dataFunctionHandle(data)
    if isa(data, 'function_handle')
        switch nargin(data)
            case 1
                fhData = @(v, t) data(v);
            case 2
                fhData = @(v, t) data(v,t);
            otherwise
                error('Incorrect number of input arguments');
        end
    else
        fhData = @(v, t) data;
    end
end

function F = stabilizedRhs(D, data)
    fhData = dataFunctionHandle(data);
    F = @(v, t, viscosity) D(v, viscosity) + fhData(v,t);
end

function F = unstabilizedRhs(D, data)
    fhData = dataFunctionHandle(data);
    F = @(v, t) D(v) + fhData(v,t);
end

function ts = rkRvTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0)
    F = stabilizedRhs(diffOpStruct.D_scheme, data);
    ts = rv.time.RungekuttaRv(F, opt.k, t0, v0, residualViscosity, diffOpStruct.D_t, opt.rkOrder);
end

function ts = rkRvBdfTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0)
    F = stabilizedRhs(diffOpStruct.D_scheme, data);
    ts = rv.time.RungekuttaRvBdf(F, opt.k, t0, v0, residualViscosity, opt.rkOrder, opt.bdfOrders);
end

function ts = rkRvMsTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0)
    F = stabilizedRhs(diffOpStruct.D_scheme, data);
    F_unstab = unstabilizedRhs(diffOpStruct.D_unstable, data);
    ts = rv.time.RungekuttaRvMultiStage(F, F_unstab, opt.k, t0, v0,...
                                        residualViscosity, diffOpStruct.D_t, opt.rkOrder);
end

function ts = rkRvMgTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0)
    F = stabilizedRhs(diffOpStruct.D_scheme, data);
    F_coarse = unstabilizedRhs(diffOpStruct.D_coarse, data);
    ts = rv.time.RungekuttaRvMultiGrid(F, F_coarse, opt.k, t0, v0,...
                                       residualViscosity, diffOpStruct.D_t, opt.rkOrder);
end

function ts = rkRvInstageTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0)
    F = stabilizedRhs(diffOpStruct.D_scheme, data);
    ts = rv.time.RungekuttaRvInstage(F, opt.k, t0, v0, residualViscosity, diffOpStruct.D_t, opt.rkOrder);
end