diff --git a/applications/solvers/combustion/plasmaReactingFoam/PhiEqn.H b/applications/solvers/combustion/plasmaReactingFoam/PhiEqn.H
index 440bb8fa..fc710b89 100644
--- a/applications/solvers/combustion/plasmaReactingFoam/PhiEqn.H
+++ b/applications/solvers/combustion/plasmaReactingFoam/PhiEqn.H
@@ -1,14 +1,22 @@
 {
+    rhoq = (rho * qc) - (eCharge * ne);
+
     solve
     (
         fvm::laplacian(Phi)
-        + fvc::Sp(rho, qc)/epsilon0
-        - fvc::Sp(eCharge, ne)/epsilon0
+        + rhoq/epsilon0
     );
 
     E = -fvc::grad(Phi);
 
     ng = p / T * (NA / R);
 
-    En = mag(E) / (ng);
+    tmp<volScalarField> tMagE (mag(E));
+    const volScalarField &magE = tMagE();
+
+    Info<< "min/max(|E|) = "
+        << min(magE).value() << ", "
+        << max(magE).value() << endl;
+
+    En = magE / (ng);
 }
diff --git a/applications/solvers/combustion/plasmaReactingFoam/createFields.H b/applications/solvers/combustion/plasmaReactingFoam/createFields.H
index 9c1e98d0..bdbde984 100644
--- a/applications/solvers/combustion/plasmaReactingFoam/createFields.H
+++ b/applications/solvers/combustion/plasmaReactingFoam/createFields.H
@@ -255,6 +255,21 @@ volScalarField ne
     mesh
 );
 
+
+Info<< "Creating field charge density\n" << endl;
+volScalarField rhoq
+(
+    IOobject
+    (
+        "rhoq",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    (rho * qc) - (eCharge * ne)
+);
+
 Info<< "Creating field gas number density\n" << endl;
 volScalarField ng("ng", p / R / T * NA);