Help with Discontinuous Galerking formulation

Hello all,
I'm facing problems with the formulation of a discontinuous galerkin.. I implemented without problemd the "CG" and now i would like to implement the "DG" and compare both... I'll put parts of my code..
My equation follow the nernst-planck diffusion equation:

$$\frac{\partial c_i}{\partial t} = \nabla \cdot (- \sum\limits_{k=1}^{Ns} D_{ik} \nabla c_i )$$

System definitions:

mesh = IntervalMesh(100,-1.0,1.0)
V = FunctionSpace(mesh.mesh, "DG", 1)
ME= MixedFunctionSpace([V] * (self._numberOfSolutes))
duList = TrialFunction(ME)
vList = TestFunctions(ME)
u = Function(ME)
uList = split(u)

Weak formulation and definition of the NonlinearVariationalProblem:

nFacet = FacetNormal(mesh)
# L = Weak statement of the equations
L = 0
## Set penalty parameter
Pe = Constant(1.0e-10)
# Cell Size
h = CellSize(mesh)
for i in range(numberOfSolutes):
    # Current unknown
    u_mid = uList[i]
    # Previous solution
    u_mid_old = u_oldList[i]
    # test function
    v = vList[i]
    LX = u_mid * v * dx - u_mid_old  * v * dx
    for k in range(numberOfSolutes):
        # Dik
        Dik = crossDiffList[k]
        # solution
        u_mid_k = uList[k]
        # composition
        q = grad(u_mid_k)
        # Weak formulation dx
        LX += dt * dot(Dik * q, grad(v)) * dx

        # as in http://www.karlin.mff.cuni.cz/~hron/fenics-tutorial/discontinuous_galerkin/doc.html
        LX += dt * Pe/avg(h) * dot(jump(u_mid_k, nFacet),jump(v, nFacet))* dS
        LX -= dt * Pe * dot(avg(grad(u_mid_k)), jump(v,nFacet)) * dS
        LX -= dt * Pe * dot(jump(u_mid_k, nFacet),avg(grad(v))) * dS
    L += LX

# Compute directional derivative about u in the direction of du (Jacobian)
a = derivative(L, u, duList)

Defining the bcs,problem and the solver:

bcs = [DirichletBC(ME.sub(i),
                                 Constant(1.0),
                                 LeftBoundary(), method="pointwise")
                     for i in range(self._numberOfSolutes)]
bcs += [DirichletBC(ME.sub(i),
                                 Constant(0.0),
                                 RightBoundary(), method="pointwise")
                     for i in range(self._numberOfSolutes)]

ffc_options = {"optimize": True, "quadrature_degree": 8}
problem = NonlinearVariationalProblem(L, u, bcs, a, form_compiler_parameters=ffc_options)
solver = NonlinearVariationalSolver(problem)
# solving the transport equations system
solver.solve()

It runs without problems, but the results are completely wrong. I reviewed many times the mathematical approach to the discontinuous galerkin but I can not find where i'm making the mistake.
I followed the instructions here , which describes a similar problem...

from here i found out that "geometric" boundary seems not to work for 1d cases.. so, i updated it to pointwise...

By reading this, i started asked myself if the problem is related to my mesh.. is DG valid for IntervalMesh at all?

Thanks in advance for any help!