How to define a given, explizit Runge-Kutta-method correctly

The implemented IMEX-Euler and IMEX-CNAB methods are working just fine, so there should somewhere be an error in the definition. In my opinion $\frac{1}{2} rhs (t_{n+1}, rhs(t_n,x_n))$ isn't implemented correctly. Thank you for your help/interest :-)
Here is my complete source code, including the IMEX-Euler and IMEX-CNAB:
from __future__ import division
# import for division, see: http://stackoverflow.    com/questions/2958684/python-division

from dolfin import *
import numpy as     np

import matplotlib as mpl
mpl.use( "agg" )
import matplotlib.pyplot as plt
from matplotlib import animation
from math import exp
from math import sin
from math import cos
# ------------------- settings, definitions
# to shut off the output
set_log_level(WARNING)

# nx = M in koto-paper
nx = 20
xmax = 1
mesh = IntervalMesh(nx, 0, xmax)

# definition of solving method
V = FunctionSpace(mesh, 'CG', 1)

T = 1
dt = 1/nx
# ------------------- boundary conditions
# Define boundary segments
class left(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-13   # tolerance for coordin    ate comparisons
        return on_boundary and abs(x[0]) <     tol

class right(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-13 
        return on_boundary and abs(x[0] - xmax) <     tol

# instances of SubDomain
Dirichlet_l = left()
Dirichlet_r = right()

# FacetFunctions, one for Neumann and one for Dirichlet boundary conditions
exterior_facet_domains_D = FacetFunction("uint", mesh)
exterior_facet_domains_D.set_all(1)

# marking of the corresponding FaceFunctions
Dirichlet_l.mark(exterior_facet_domains_D, 2)
Dirichlet_r.mark(exterior_facet_domains_D, 3)

# definition of g_N and g_D, corresponding Functions
g_0 = Expression("exp(t) * sin(t)", t=0)
g_1 = Expression("exp(t-1)*sin(t+1)", t=0)

# definition of Dirichlet boundary conditions, Neumann     boundaries in variational form
bcs = [DirichletBC(V, g_0, exterior_facet_domains_D, 2),
       DirichletBC(V, g_1, exterior_facet_domains_D, 3)]
# ------------------- boundary conditions

# ------------------- RHS, exact solution
# RHS
def rhs(u):
  return u - u*u + f

def rhs_old(u):
  return u - u*u + f_old

def rhs_old_2(u):
  return u - u*u + f_old_2

f = Expression("3*exp(t-x[0])*cos(t+x[0])+exp(2*(t-x[0]))*(pow(sin(t+x[0]),2))", t=0)
f_old = Expression("3*exp(t-x[0])*cos(t+x[0])+exp(2*(t-x[0]))*(pow(sin(t+x[0]),2))", t=0)
f_old_2 = Expression("3*exp(t-x[0])*cos(t+x[0])+exp(2*(t-x[0]))*(pow(sin(t+x[0]),2))", t=0)

u_exact = Expression("exp(t-x[0])*sin(t+x[0])",t=0)

# initial values
u_start = Expression("exp(-x[0])*sin(x[0])")
# ------------------- RHS, exact solution

# ------------------- IMEX
# ---------- Euler
m1_old = interpolate(u_start,V)
m1 = TrialFunction(V)
n1 = TestFunction(V)

B1 = inner(m1,n1)*dx + dt*inner(nabla_grad(m1), nabla_grad(n1))*dx
G1 = inner(dt*rhs_old(m1_old) + m1_old, n1)*dx

m1 = Function(V)
m1 = interpolate(u_start,V)
# ---------- CNAB
m2_old = interpolate(u_start,V)
m2_old_2 = interpolate(u_start,V)
m2 = TrialFunction(V)
n2 = TestFunction(V)

B2 = inner(m2,n2)*dx + dt*1/2*inner(nabla_grad(m2), nabla_grad(n2))*dx
G2 = - dt*1/2*inner(nabla_grad(m2_old), nabla_grad(n2))*dx + inner(dt*3/2*rhs_old(m2_old) - dt*1/2* rhs_old_2(m2_old_2) + m2_old, n2)*dx

m2 = Function(V)
m2 = interpolate(u_start,V)
# ---------- trapezoidal, tpz
m3_old = interpolate(u_start,V)
m3 = TrialFunction(V)
n3 = TestFunction(V)

k3_1 = rhs_old(m3_old)
k3_2 = rhs(k3_1)

B3 = inner(m3,n3)*dx + dt*1/2*inner(nabla_grad(m3), nabla_grad(n3))*dx
G3 = - dt*1/2*inner(nabla_grad(m3_old), nabla_grad(n3))*dx + inner(dt*1/2*(k3_1 + k3_2) + m3_old, n3)*dx

m3 = Function(V)
m3 = interpolate(u_start,V)
# ------------------- IMEX

# ------------------- 
norm_IMEX = 0
norm_CNAB = 0
norm_tpz = 0

# as first step is just interpolating, it can be saved directly
u_exact_int = interpolate(u_exact,V)

# ---------- update of time dependent data
t = dt
g_0.t = t
g_1.t = t
f.t = t
f_old.t = t-dt
f_old_2.t = t-2*dt
u_exact.t=t
u_exact_int = Function(V)

# ---------- time steps
while t <= T+dt:

    # ---------- IMEX
    # Euler
    solve(B1 == G1, m1, bcs = bcs)
    #tpz
    if(t >= 2*dt):
      solve(B3 == G3, m3, bcs = bcs)
    else:
      m3.assign(m1)
    # CNAB
    if(t >= 2*dt):
      solve(B2 == G2, m2, bcs = bcs)
    else:
      m2.assign(m1)
    # ---------- computation of error     
    error = errornorm(u_exact, m3,'l2',5)
    if(error>norm_tpz):
       norm_tpz = error

    error = errornorm(u_exact, m1,'l2',5)
    if(error>norm_IMEX ):
       norm_IMEX  = error

    error = errornorm(u_exact, m2,'l2',5)
    if(error>norm_CNAB):
       norm_CNAB = error

    # ---------- interpolation of exact solution
    u_exact_int = interpolate(u_exact,V)

    # ---------- assigning of solution from previous time step
    m1_old.assign(m1)
    m2_old_2.assign(m2_old)
    m2_old.assign(m2)
    m3_old.assign(m3)

    # ---------- update of time dependent data
    t += dt
    g_0.t = t
    g_1.t = t
    f.t = t
    f_old.t = t-dt
    f_old_2.t = t-2*dt
    u_exact.t=t

#  ------------------- ploting, output of errors
print 'L2-Error IMEX        = %.2e' % (norm_IMEX)
print 'L2-Error IMEX CNAB   = %.2e' % (norm_CNAB)
print 'L2-Error IMEX tpz    =%.12f' % (norm_tpz)