def METHOD_NAME(timeintegrator_class, timeintegrator_options, refinement=1):
    """
    Run test for the given time integrator
    """
    mesh = UnitSquareMesh(2, 2)
    # TODO create ThetisMesh and push the bnd def there
    bnd_len = compute_boundary_length(mesh)
    mesh.boundary_len = bnd_len
    p1 = get_functionspace(mesh, 'CG', 1)
    fs = MixedFunctionSpace([p1, p1])
    alpha = 2*numpy.pi
    mode = 'explicit'
    cfl_coeff = timeintegrator_class.cfl_coeff if hasattr(timeintegrator_class, 'cfl_coeff') else None
    if (cfl_coeff == CFL_UNCONDITIONALLY_STABLE
            or isinstance(timeintegrator_class, DIRKGeneric)):
        mode = 'implicit'
    if (IMEXGeneric in timeintegrator_class.__bases__):
        mode = 'imex'
    equation = SimpleODEEquation(fs, alpha, mode=mode)
    solution = Function(fs, name='solution')
    solution.sub(0).assign(Constant(0))
    solution.sub(1).assign(Constant(1))
    fields = {}
    end_time = 1.0
    base_dt = 0.01
    ntimesteps = int(numpy.round(end_time/base_dt*refinement))
    dt = end_time/ntimesteps
    print('Running refinement {:2d}, dt = {:.6f}'.format(refinement, dt))
    times = numpy.zeros((ntimesteps+1, ))
    values = numpy.zeros((ntimesteps+1, 2))
    ti = timeintegrator_class(equation, solution, fields, dt, timeintegrator_options(), {})
    ti.initialize(solution)
    simulation_time = 0
    sol_a, sol_b = solution.split()
    values[0, :] = sol_a.dat.data[0], sol_b.dat.data[0]
    for i in range(ntimesteps):
        simulation_time = (i+1)*dt
        ti.advance(simulation_time)
        times[i+1] = simulation_time
        values[i+1, :] = sol_a.dat.data[0], sol_b.dat.data[0]
    assert abs(times[-1] - end_time) < 1e-16
    exact_sol = numpy.vstack((numpy.sin(alpha*times), numpy.cos(alpha*times))).T
    l2_err = numpy.sqrt(numpy.mean((values - exact_sol)**2))
    return l2_err