def METHOD_NAME(M, G, 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())
    while (t < params["finish_time"]):
        t += dt
        step += 1
        rhs = action(A_r, state)
        # Solve the shallow water equations.
        solve(A == rhs, tmpstate, annotate=annotate)
        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
    return state  # return the state at the final time