diff --git a/applications/solvers/combustion/plasmaReactingFoam/PhiEqn.H b/applications/solvers/combustion/plasmaReactingFoam/PhiEqn.H
index 8b78a8b2..dc857322 100644
--- a/applications/solvers/combustion/plasmaReactingFoam/PhiEqn.H
+++ b/applications/solvers/combustion/plasmaReactingFoam/PhiEqn.H
@@ -7,4 +7,5 @@
     );
 
     E = -fvc::grad(Phi);
+    En = mag(E) / (ng);
 }
diff --git a/applications/solvers/combustion/plasmaReactingFoam/YEqn.H b/applications/solvers/combustion/plasmaReactingFoam/YEqn.H
index bc9726d0..4e2177be 100644
--- a/applications/solvers/combustion/plasmaReactingFoam/YEqn.H
+++ b/applications/solvers/combustion/plasmaReactingFoam/YEqn.H
@@ -15,18 +15,38 @@ tmp<fv::convectionScheme<scalar> > mvConvection
     label inertIndex = -1;
     volScalarField Yt(0.0*Y[0]);
 
+    forAll(rho, celli)
+    {
+        const scalar Eni = En[celli];
+        De[celli] = DeN.value(Eni);
+        mue[celli] = mueN.value(Eni);
+        Te[celli] = TeOfEn.value(Eni);
+    }
+    De.correctBoundaryConditions();
+    mue.correctBoundaryConditions();
+    Te.correctBoundaryConditions();
+
     forAll(Y, i)
     {
         if (Y[i].name() == electronSpecie)
         {
+            volScalarField& Yi = Y[i];
+
+            ve = - linearInterpolate(mue*E/ng) & mesh.Sf();
+
+            tmp<fvScalarMatrix> electronR(
+                    new fvScalarMatrix(ne, dimless/dimTime));
+            electronR->source() = reaction->R(Yi)->source();
+
             fvScalarMatrix neEqn
                 (
                  fvm::ddt(ne)
                  + mvConvection->fvmDiv(ve, ne)
-                 - mvConvection->fvmDiv(fvc::interpolate(De/Te*fvc::grad(Te)) & mesh.Sf(), ne)
-                 - fvm::laplacian(De, ne)
+                 - mvConvection->fvmDiv(fvc::interpolate(De/ng/Te*fvc::grad(Te)) & mesh.Sf(), ne)
+                 - fvm::laplacian(De/ng, ne)
                  ==
-                 fvOptions(ne)
+                 electronR
+                 + fvOptions(ne)
                 );
 
             neEqn.relax();
diff --git a/applications/solvers/combustion/plasmaReactingFoam/createFields.H b/applications/solvers/combustion/plasmaReactingFoam/createFields.H
index 185425ab..06cd1407 100644
--- a/applications/solvers/combustion/plasmaReactingFoam/createFields.H
+++ b/applications/solvers/combustion/plasmaReactingFoam/createFields.H
@@ -22,6 +22,9 @@ dimensionedScalar k
     physicalProperties.lookup("k")
 );
 
+CSV<scalar> mueN("mueN", physicalProperties);
+CSV<scalar> DeN("DeN", physicalProperties);
+CSV<scalar> TeOfEn("TeOfEn", physicalProperties);
 
 Info<< "Reading field Phi\n" << endl;
 volScalarField Phi
@@ -172,7 +175,7 @@ Info<< "Creating field gas number density\n" << endl;
 volScalarField ng("ng", p / R / T * NA);
 
 Info<< "Creating field reduced electric field\n" << endl;
-volScalarField En ("En", mag(E) / (ne+ng));
+volScalarField En ("En", mag(E) / (ng));
 
 Info<< "Creating field electron mobility\n" << endl;
 volScalarField mue
@@ -216,6 +219,8 @@ volScalarField Te
     mesh
 );
 
+
+/*
 BE_IX::Bolos bolos;
 
 bolos.readCollisions("./LXCat-June2013.txt");
@@ -231,6 +236,7 @@ bolos.presolve();
     De.correctBoundaryConditions();
     mue.correctBoundaryConditions();
 }
+*/
 
 Info<< "Calculating face flux field ve\n" << endl;
 surfaceScalarField ve
@@ -243,6 +249,6 @@ surfaceScalarField ve
         IOobject::NO_READ,
         IOobject::NO_WRITE
     ),
-    linearInterpolate(mue*E) & mesh.Sf()
+    - linearInterpolate(mue*E/ng) & mesh.Sf()
 );
 
diff --git a/applications/solvers/combustion/plasmaReactingFoam/plasmaReactingFoam.C b/applications/solvers/combustion/plasmaReactingFoam/plasmaReactingFoam.C
index 6b2b3950..bdc28195 100644
--- a/applications/solvers/combustion/plasmaReactingFoam/plasmaReactingFoam.C
+++ b/applications/solvers/combustion/plasmaReactingFoam/plasmaReactingFoam.C
@@ -36,6 +36,8 @@ Description
 #include "pimpleControl.H"
 #include "fvIOoptionList.H"
 
+#include "CSV.H"
+
 #include "bolos.h"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -69,7 +71,6 @@ int main(int argc, char *argv[])
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
         #include "PhiEqn.H"
-        En = mag(E) / (ne+ng);
         #include "rhoEqn.H"
 
         while (pimple.loop())