def METHOD_NAME(M, G, rhs_contr, ufl, ufr, state, params, annotate=True):
    '''Solve M*dstate/dt = G*state using a theta scheme.'''
    A = M + params["theta"] * params["dt"] * G
    A_r = M - (1 - params["theta"]) * params["dt"] * G
    u_out, p_out = output_files(params["basename"])
    M_u_out, v_out, u_out_state = u_output_projector(state.function_space())
    M_p_out, q_out, p_out_state = p_output_projector(state.function_space())
    # Project the solution to P1 for visualisation.
    rhs = assemble(inner(v_out, state.split()[0]) * dx)
    solve(M_u_out, u_out_state.vector(), rhs, "cg", "sor", annotate=False)
    # Project the solution to P1 for visualisation.
    rhs = assemble(inner(q_out, state.split()[1]) * dx)
    solve(M_p_out, p_out_state.vector(), rhs, "cg", "sor", annotate=False)
    u_out << u_out_state
    p_out << p_out_state
    t = params["current_time"]
    dt = params["dt"]
    step = 0
    tmpstate = Function(state.function_space())
    j = 0
    (u_j, p_j) = split(state)
    j += 0.5 * dt * assemble(dot(u_j, u_j) * dx)
    while (t < params["finish_time"]):
        t += dt
        ufl.t = t - (1.0 - params["theta"]) * dt  # Update time for the Boundary condition expression
        step += 1
        rhs = action(A_r, state) + params["dt"] * rhs_contr
        # Solve the shallow water equations.
        solve(A == rhs, tmpstate, annotate=annotate)
        state.assign(tmpstate, annotate=annotate)
        # We solve twice, to make sure that Functional handles multiple-solves-per-timestep correctly
        solve(A == rhs, tmpstate, annotate=annotate)
        # solve(A, state.vector(), rhs, "preonly", "lu")
        state.assign(tmpstate, annotate=annotate)
        if step % params["dump_period"] == 0:
            # Project the solution to P1 for visualisation.
            rhs = assemble(inner(v_out, state.split()[0]) * dx)
            solve(M_u_out, u_out_state.vector(), rhs, "cg", "sor", annotate=False)
            # Project the solution to P1 for visualisation.
            rhs = assemble(inner(q_out, state.split()[1]) * dx)
            solve(M_p_out, p_out_state.vector(), rhs, "cg", "sor", annotate=False)
            u_out << u_out_state
            p_out << p_out_state
        if t >= params["finish_time"]:
            quad_weight = 0.5
        else:
            quad_weight = 1.0
        (u_j, p_j) = split(state)
        j += quad_weight * dt * assemble(dot(u_j, u_j) * dx)
    return j, state  # return the state and the functional's time integral contribution at the final time