diff --git a/.gitignore b/.gitignore
index 0a9a292..ca181ee 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,4 +1,4 @@
-DSMC_Neutrals
+PPartiC
 mod/
 obj/
 doc/user_manual/
diff --git a/makefile b/makefile
index b10a4ac..3dfb28d 100644
--- a/makefile
+++ b/makefile
@@ -12,20 +12,24 @@ JSONINC := $(JSONDIR)/include
 INCDIR := 
 INCDIR += $(JSONINC)
 
+VAR := ""
+
 # compile flags
-FCFLAGS := -fopenmp -Ofast -J $(MODDIR) -I $(INCDIR) -Wall -fno-stack-arrays -g
+FCFLAGS := -fopenmp -Ofast -J $(MODDIR) -I $(INCDIR) -Wall -g
 
 #Output file
-OUTPUT = DSMC_Neutrals
+OUTPUT = PPartiC
 
 export
 
-$(OUTPUT): src.o 
+all: $(OUTPUT)
+
+$(OUTPUT):  src.o 
 
 src.o:
 	@mkdir -p $(MODDIR)
 	@mkdir -p $(OBJDIR)
-	$(MAKE) -C src all
+	$(MAKE) -C src $(OUTPUT)
 
 clean:
 	rm -f $(OUTPUT)
diff --git a/runs/cylFlow/input.json b/runs/cylFlow/input.json
index 04814f7..440fb58 100644
--- a/runs/cylFlow/input.json
+++ b/runs/cylFlow/input.json
@@ -24,14 +24,15 @@
     {"name": "Nozzle", "species": "Argon", "flow": 10.0, "units": "sccm", "v": 300.0, "T": [300.0, 300.0, 300.0], "n": [1, 0, 0], "physicalSurface": 1}
   ],
   "reference": {
-	"density": 1.0e20,
-	"mass": 6.633e-26,
-	"temperature": 300.0,
-	"radius": 1.88e-10
+    "density": 1.0e20,
+    "mass": 6.633e-26,
+    "temperature": 300.0,
+    "radius": 1.88e-10
   },
   "case": {
-    "tau":  6.e-7,
-    "time": 3.e-3
+    "tau":  6.0e-7,
+    "time": 3.0e-3,
+    "solver": "2DCylNeutral"
   },
   "interactions": {
     "folderCollisions": "./data/collisions/",
diff --git a/src/DSMC_Neutrals.f95 b/src/PPartiC.f95
similarity index 77%
rename from src/DSMC_Neutrals.f95
rename to src/PPartiC.f95
index f71bdbb..e7db5c7 100644
--- a/src/DSMC_Neutrals.f95
+++ b/src/PPartiC.f95
@@ -1,5 +1,5 @@
-!  PPartiC neutral DSMC main program.
-PROGRAM DSMC_Neutrals
+!  PPartiC neutral PIC main program.
+PROGRAM PPartiC
   USE moduleInput
   USE moduleErrors
   USE OMP_LIB
@@ -10,28 +10,32 @@ PROGRAM DSMC_Neutrals
   USE moduleMesh
   IMPLICIT NONE
 
-  ! t = time step; n = number of particle; e = number of element; i = number of inject
+  ! t = time step
   INTEGER:: t = 0
+  ! arg1 = Input argument 1 (input file)
   CHARACTER(200):: arg1
+  ! inputFile = path+name of input file
   CHARACTER(:), ALLOCATABLE:: inputFile
 
-
   tStep = omp_get_wtime()
   !Gets the input file
   CALL getarg(1, arg1)
   inputFile = TRIM(arg1)
   !If no input file is declared, a critical error is called
-  IF (inputFile == "") CALL criticalError("No input file provided", "DSMC_Neutrals")
+  IF (inputFile == "") CALL criticalError("No input file provided", "PPartiC")
 
   !Reads the json configuration file
   CALL readConfig(inputFile)
 
+  CALL verboseError("Calculating initial EM field...")
+  CALL doEMField()
+
   tStep = omp_get_wtime() - tStep
 
   !Output initial state
   CALL doOutput(t)
   CALL verboseError('Starting main loop...')
-  !$OMP PARALLEL PRIVATE(t) DEFAULT(SHARED) 
+  !$OMP PARALLEL DEFAULT(SHARED) 
   DO t = 1, tmax
     !Insert new particles and push them
     !$OMP SINGLE
@@ -46,18 +50,18 @@ PROGRAM DSMC_Neutrals
 
     !$OMP SINGLE
     tPush = omp_get_wtime() - tPush
+    !Collisions
     tColl = omp_get_wtime()
     !$OMP END SINGLE
 
-    !Collisions
     CALL doCollisions()
 
     !$OMP SINGLE
     tColl = omp_get_wtime() - tColl
+    !Reset particles
     tReset = omp_get_wtime()
     !$OMP END SINGLE
 
-    !Reset particles
     CALL doReset()
 
     !$OMP SINGLE
@@ -71,6 +75,13 @@ PROGRAM DSMC_Neutrals
 
     !$OMP SINGLE
     tWeight = omp_get_wtime() - tWeight
+    tEMField = omp_get_wtime()
+    !$OMP END SINGLE
+
+    CALL doEMField()
+
+    !$OMP SINGLE
+    tEMField = omp_get_wtime() - tEMField
     tStep = omp_get_wtime() - tStep
     !Output data
     CALL doOutput(t)
@@ -79,5 +90,5 @@ PROGRAM DSMC_Neutrals
   END DO
   !$OMP END PARALLEL
 
-END PROGRAM DSMC_Neutrals
+END PROGRAM PPartiC
 
diff --git a/src/makefile b/src/makefile
index b0dfff0..f9b4928 100644
--- a/src/makefile
+++ b/src/makefile
@@ -1,17 +1,16 @@
-OBJECTS = $(OBJDIR)/$(OUTPUT).o \
-	        $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleCompTime.o $(OBJDIR)/moduleSolver.o \
+OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleCompTime.o $(OBJDIR)/moduleSolver.o \
 			  	$(OBJDIR)/moduleSpecies.o $(OBJDIR)/moduleInject.o $(OBJDIR)/moduleInput.o \
 			  	$(OBJDIR)/moduleErrors.o $(OBJDIR)/moduleList.o $(OBJDIR)/moduleOutput.o \
-					$(OBJDIR)/moduleMeshCyl.o	$(OBJDIR)/moduleMeshCylRead.o	$(OBJDIR)/moduleMeshCylBoundary.o \
-					$(OBJDIR)/moduleBoundary.o $(OBJDIR)/moduleCaseParam.o $(OBJDIR)/moduleRefParam.o \
-					$(OBJDIR)/moduleCollisions.o $(OBJDIR)/moduleTable.o $(OBJDIR)/moduleParallel.o
-
+			  	$(OBJDIR)/moduleMeshCyl.o	$(OBJDIR)/moduleMeshCylRead.o	$(OBJDIR)/moduleMeshCylBoundary.o \
+			  	$(OBJDIR)/moduleBoundary.o $(OBJDIR)/moduleCaseParam.o $(OBJDIR)/moduleRefParam.o \
+			  	$(OBJDIR)/moduleCollisions.o $(OBJDIR)/moduleTable.o $(OBJDIR)/moduleParallel.o \
+			  	$(OBJDIR)/moduleEM.o
 
 all: $(OUTPUT)
-	$(FC) $(FCFLAGS) -o $(TOPDIR)/$(OUTPUT) $(OBJECTS) $(JSONLIB)
 
 $(OUTPUT): modules.o $(OUTPUT).f95
 	$(FC) $(FCFLAGS) -o $(OBJDIR)/$(OUTPUT).o -c $(OUTPUT).f95
+	$(FC) $(FCFLAGS) -o $(TOPDIR)/$(OUTPUT) $(OBJECTS) $(OBJDIR)/$(OUTPUT).o $(JSONLIB) -L/usr/local/lib -llapack -lopenblas
 
 modules.o:
 	$(MAKE) -C modules all
diff --git a/src/modules/makefile b/src/modules/makefile
index 0cf8c24..e091015 100644
--- a/src/modules/makefile
+++ b/src/modules/makefile
@@ -1,9 +1,9 @@
 
 OBJS = moduleCaseParam.o moduleCompTime.o moduleList.o\
-			 moduleMesh.o moduleMeshCyl.o moduleMeshCylBoundary.o\
-			 moduleMeshCylRead.o moduleOutput.o moduleInput.o \
-			 moduleSolver.o moduleCollisions.o moduleTable.o \
-			 moduleParallel.o
+		   moduleMesh.o moduleMeshCyl.o moduleMeshCylBoundary.o\
+		   moduleMeshCylRead.o moduleOutput.o moduleInput.o \
+		   moduleSolver.o moduleCollisions.o moduleTable.o \
+		   moduleParallel.o moduleEM.o
 
 all: $(OBJS)
 
@@ -40,11 +40,14 @@ moduleRefParam.o: moduleConstParam.o moduleRefParam.f95
 moduleSpecies.o: moduleCaseParam.o moduleSpecies.f95
 	$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
 
-moduleSolver.o: moduleSpecies.o moduleRefParam.o moduleMesh.o moduleOutput.o moduleSolver.f95
+moduleSolver.o: moduleEM.o moduleSpecies.o moduleRefParam.o moduleMesh.o moduleOutput.o moduleSolver.f95
 	$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
 
 moduleTable.o: moduleErrors.o moduleTable.f95
 	$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
 
+moduleEM.o: moduleSpecies.o moduleMesh.o moduleEM.f95
+	$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
+
 %.o: %.f95
 	$(FC) $(FCFLAGS) -c $< -o $(OBJDIR)/$@
diff --git a/src/modules/moduleCompTime.f95 b/src/modules/moduleCompTime.f95
index a45a07d..28273ce 100644
--- a/src/modules/moduleCompTime.f95
+++ b/src/modules/moduleCompTime.f95
@@ -5,6 +5,7 @@ MODULE moduleCompTime
   REAL(8):: tReset=0.D0
   REAL(8):: tColl=0.D0
   REAL(8):: tWeight=0.D0
+  REAL(8):: tEMField=0.D0
 
 END MODULE moduleCompTime
 
diff --git a/src/modules/moduleConstParam.f95 b/src/modules/moduleConstParam.f95
index a335a53..634f845 100644
--- a/src/modules/moduleConstParam.f95
+++ b/src/modules/moduleConstParam.f95
@@ -6,9 +6,10 @@ MODULE moduleConstParam
 
   REAL(8), PARAMETER:: PI = 4.D0*DATAN(1.D0) !number pi
   REAL(8), PARAMETER:: sccm2atomPerS = 4.5D17 !sccm to atom s^-1
-  REAL(8), PARAMETER:: eV2J = 1.60217662D-19 !Electron volt to Joule (elementary charge)
+  REAL(8), PARAMETER:: qe = 1.60217662D-19 !Elementary charge
   REAL(8), PARAMETER:: kb = 1.38064852D-23 !Boltzmann constants SI
-  ! REAL(8), PARAMETER:: eps_0 = 8.8542D-12 !Epsilon_0 (Not used in Neutrals)
+  REAL(8), PARAMETER:: eV2J = qe !Electron volt to Joule conversion
+  REAL(8), PARAMETER:: eps_0 = 8.8542D-12 !Epsilon_0
 
 END MODULE moduleConstParam
 
diff --git a/src/modules/moduleEM.f95 b/src/modules/moduleEM.f95
new file mode 100644
index 0000000..996bd33
--- /dev/null
+++ b/src/modules/moduleEM.f95
@@ -0,0 +1,133 @@
+!Module to solve the electromagnetic field
+MODULE moduleEM
+  IMPLICIT NONE
+
+  TYPE:: boundaryEM
+    CHARACTER(:), ALLOCATABLE:: typeEM
+    INTEGER:: physicalSurface
+    REAL(8):: potential
+
+    CONTAINS
+      PROCEDURE, PASS:: apply
+
+  END TYPE boundaryEM
+
+  INTEGER:: nBoundaryEM
+  TYPE(boundaryEM), ALLOCATABLE:: boundEM(:)
+
+  !Information of charge and reference parameters for rho vector
+  REAL(8), ALLOCATABLE:: qSpecies(:)
+
+  CONTAINS
+    SUBROUTINE apply(self, edge)
+      USE moduleMesh
+      IMPLICIT NONE
+
+      CLASS(boundaryEM), INTENT(in):: self
+      CLASS(meshEdge):: edge
+      INTEGER, ALLOCATABLE:: nodes(:)
+      INTEGER:: n
+
+      nodes = edge%getNodes()
+
+      DO n = 1, SIZE(nodes)
+        SELECT CASE(self%typeEM)
+        CASE ("dirichlet")
+          mesh%K(nodes(n), :)        = 0.D0
+          mesh%K(nodes(n), nodes(n)) = 1.D0
+        
+          mesh%nodes(nodes(n))%obj%emData%type = self%typeEM
+          mesh%nodes(nodes(n))%obj%emData%phi  = self%potential
+
+        END SELECT
+
+      END DO
+
+    END SUBROUTINE
+
+    PURE FUNCTION gatherElecField(part) RESULT(elField)
+      USE moduleSpecies
+      USE moduleMesh
+      IMPLICIT NONE
+
+      TYPE(particle), INTENT(in):: part
+      REAl(8):: xi(1:3) !Logical coordinates of particle in element
+      REAL(8):: elField(1:3)
+
+      elField = 0.D0
+
+      xi  = part%xLog
+
+      elField = mesh%vols(part%e_p)%obj%gatherEF(xi)
+
+    END FUNCTION gatherElecField
+
+    SUBROUTINE assembleSourceVector(vectorF)
+      USE moduleMesh
+      USE moduleRefParam
+      IMPLICIT NONE
+
+      REAL(8), INTENT(out):: vectorF(1:mesh%numNodes)
+      REAL(8), ALLOCATABLE:: localF(:)
+      INTEGER, ALLOCATABLE:: nodes(:)
+      REAL(8), ALLOCATABLE:: rho(:)
+      INTEGER:: nNodes
+      INTEGER:: e, i, ni
+      CLASS(meshNode), POINTER:: node
+
+      !$OMP SINGLE
+      vectorF = 0.D0
+      !$OMP END SINGLE
+
+      !$OMP SINGLE
+      DO e = 1, mesh%numVols
+        nodes  = mesh%vols(e)%obj%getNodes()
+        nNodes = SIZE(nodes)
+        !Calculates charge density (rho) in element nodes
+        ALLOCATE(rho(1:nNodes))
+        rho = 0.D0
+        DO i = 1, nNodes
+          ni = nodes(i)
+          node => mesh%nodes(ni)%obj
+          rho(i) = DOT_PRODUCT(qSpecies(:), node%output(:)%den/(vol_ref*node%v*n_ref))
+
+        END DO
+
+        !If there is charge in the nodes, compute the local F vector
+        IF (ANY(rho > 0.D0)) THEN
+          !Calculates local F vector
+          localF = mesh%vols(e)%obj%elemF(rho)
+
+          !Assign local F to global F
+          DO i = 1, nNodes
+            ni = nodes(i)
+            vectorF(ni) = vectorF(ni) + localF(i)
+
+          END DO
+
+          DEALLOCATE(localF)
+
+        END IF
+
+        DEALLOCATE(nodes, rho)
+
+      END DO
+      !$OMP END SINGLE
+
+      !Apply boundary conditions
+      !$OMP DO
+      DO i = 1, mesh%numNodes
+        node => mesh%nodes(i)%obj
+
+        SELECT CASE(node%emData%type)
+        CASE ("dirichlet")
+          vectorF(i) = node%emData%phi
+
+        END SELECT
+
+      END DO
+      !$OMP END DO
+
+    END SUBROUTINE assembleSourceVector
+
+END MODULE moduleEM
diff --git a/src/modules/moduleErrors.f95 b/src/modules/moduleErrors.f95
index 1d61142..0a95db8 100644
--- a/src/modules/moduleErrors.f95
+++ b/src/modules/moduleErrors.f95
@@ -37,7 +37,6 @@ MODULE moduleErrors
       errorMsg  = msg
 
       WRITE (*, '(A)') errorMsg
-      WRITE (*, *)
 
     END SUBROUTINE verboseError
 
diff --git a/src/modules/moduleInject.f95 b/src/modules/moduleInject.f95
index 9b87133..2fe2e6f 100644
--- a/src/modules/moduleInject.f95
+++ b/src/modules/moduleInject.f95
@@ -25,24 +25,14 @@ MODULE moduleInject
   TYPE(injectGeneric), ALLOCATABLE:: inject(:)
 
   CONTAINS
-    !Injection of particles
-    SUBROUTINE doInjects()
-      IMPLICIT NONE
-
-      INTEGER:: i
-
-      DO i=1, nInject
-        CALL inject(i)%addParticles()
-      END DO
-
-    END SUBROUTINE doInjects
-
-    SUBROUTINE initInject(self, i, v, n, T, flow, pt, physicalSurface)
+    !Initialize an injection of particles
+    SUBROUTINE initInject(self, i, v, n, T, flow, units, pt, physicalSurface)
       USE moduleMesh
       USE moduleMeshCyl
       USE moduleRefParam
       USE moduleConstParam
       USE moduleSpecies
+      USE moduleErrors
       IMPLICIT NONE
 
       CLASS(injectGeneric), INTENT(inout):: self
@@ -50,6 +40,7 @@ MODULE moduleInject
       REAL(8), INTENT(in):: v, n(1:3), T(1:3)
       INTEGER, INTENT(in):: pt, physicalSurface
       REAL(8), INTENT(in):: flow
+      CHARACTER(:), ALLOCATABLE, INTENT(in):: units
       INTEGER:: e, et
       INTEGER:: phSurface(1:mesh%numEdges)
 
@@ -57,8 +48,21 @@ MODULE moduleInject
       self%vMod       = v/v_ref
       self%n          = n
       self%T          = T/T_ref
-      self%nParticles = INT(flow*sccm2atomPerS*tau*ti_ref/species(pt)%obj%weight)
-      self%pt         = pt
+      SELECT CASE(units)
+      CASE ("sccm")
+        !Standard cubic centimeter per minute
+        self%nParticles = INT(flow*sccm2atomPerS*tau*ti_ref/species(pt)%obj%weight)
+
+      CASE ("A")
+        !Input current in Ampers
+        self%nParticles = INT(flow*tau*ti_ref/(qe*species(pt)%obj%weight))
+
+      CASE DEFAULT
+        CALL criticalError("No support for units: " // units, 'initInject')
+        
+      END SELECT
+      IF (self%nParticles == 0) CALL criticalError("The number of particles for inject is 0.", 'initInject')
+      self%pt  = pt
 
       self%v   = self%v * self%n
       self%vTh = DSQRT(self%T/species(self%pt)%obj%m)
@@ -92,6 +96,19 @@ MODULE moduleInject
 
     END SUBROUTINE
 
+    !Injection of particles
+    SUBROUTINE doInjects()
+    
+      IMPLICIT NONE
+
+      INTEGER:: i
+
+      DO i=1, nInject
+        CALL inject(i)%addParticles()
+      END DO
+
+    END SUBROUTINE doInjects
+
     !Random velocity from maxwellian distribution
     REAL(8) FUNCTION vBC(u, vth)
       USE moduleConstParam, ONLY: PI
@@ -119,18 +136,28 @@ MODULE moduleInject
       CLASS(injectGeneric), INTENT(in):: self
       REAL(8):: randomX
       INTEGER:: j
-      REAL(8):: randomPos
-      INTEGER:: nMin, nMax !Min and Max index in partInj array
+      INTEGER, SAVE:: nMin, nMax !Min and Max index in partInj array
       INTEGER:: n
       CLASS(meshEdge), POINTER:: randomEdge
 
       !Insert particles
+      !$OMP SINGLE
       nMin = SUM(inject(1:(self%id-1))%nParticles) + 1
-      nMax = nMin + self%nParticles
+      nMax = nMin + self%nParticles - 1
       !Assign particle type
       partInj(nMin:nMax)%pt = self%pt
       !Assign weight to particle.
       partInj(nMin:nMax)%weight = species(self%pt)%obj%weight
+      !Particle is considered to be outside the domain
+      partInj(nMin:nMax)%n_in = .FALSE.
+      !Assign charge/mass to charged particle.
+      SELECT TYPE(sp => species(self%pt)%obj)
+      TYPE IS(speciesCharged)
+        partInj(nMin:nMax)%qm = sp%qm
+
+      END SELECT
+      !$OMP END SINGLE
+
       !$OMP DO
       DO n = nMin, nMax
         !Select edge randomly from which inject particle
@@ -143,8 +170,7 @@ MODULE moduleInject
 
         randomEdge => mesh%edges(self%edges(j))%obj
         !Random position in edge
-        randomPos = RAND()
-        partInj(n)%r = randomEdge%randPos(randomPos)
+        partInj(n)%r = randomEdge%randPos()
         !Volume associated to the edge:
         IF (DOT_PRODUCT(self%n, randomEdge%normal) >= 0.D0) THEN
           partInj(n)%e_p = randomEdge%e1%n
@@ -159,7 +185,7 @@ MODULE moduleInject
                           vBC(self%v(3), self%vTh(3)) /)
 
         !Push new particle
-        CALL push(partInj(n))
+        CALL solver%push(partInj(n))
         !Assign cell to new particle
         CALL mesh%vols(partInj(n)%e_p)%obj%findCell(partInj(n))
 
diff --git a/src/modules/moduleInput.f95 b/src/modules/moduleInput.f95
index 926147b..e489a75 100644
--- a/src/modules/moduleInput.f95
+++ b/src/modules/moduleInput.f95
@@ -44,6 +44,9 @@ MODULE moduleInput
       CALL verboseError('Reading Geometry...')
       CALL readGeometry(config)
 
+      !Read boundary for EM field
+      CALL readEMBoundary(config)
+
       !Read injection of particles
       CALL verboseError('Reading Interactions between species...')
       CALL readInject(config)
@@ -63,7 +66,7 @@ MODULE moduleInput
       IMPLICIT NONE
 
       TYPE(json_file), INTENT(inout):: config
-      LOGICAL:: found
+      LOGICAL:: found, found_r
       CHARACTER(:), ALLOCATABLE:: object
 
       object = 'reference'
@@ -77,15 +80,25 @@ MODULE moduleInput
       CALL config%get(object // '.temperature', T_ref, found)
       IF (.NOT. found) CALL criticalError('Reference temperature not found','readReference')
 
-      CALL config%get(object // '.radius', r_ref, found)
-      IF (.NOT. found) CALL criticalError('Reference radius not found','readReference')
-
+      CALL config%get(object // '.radius', r_ref, found_r)
+      
       !Derived parameters
       v_ref     = DSQRT(kb*T_ref/m_ref)  !reference velocity
-      sigma_ref = PI*(r_ref+r_ref)**2    !reference cross section
-      L_ref     = 1.D0/(sigma_ref*n_ref) !mean free path
+      !TODO: Make this solver dependent
+      IF (found_r) THEN
+        L_ref = 1.D0/(sigma_ref*n_ref)     !mean free path
+        sigma_ref = PI*(r_ref+r_ref)**2    !reference cross section
+
+      ELSE
+        L_ref = DSQRT(kb*T_ref*eps_0/n_ref)/qe !Debye length
+        !TODO: Obtain right sigma_ref for PIC case
+        sigma_ref = PI*(4.D-10)**2        !reference cross section
+
+      END IF
       ti_ref    = L_ref/v_ref            !reference time
       Vol_ref   = L_ref**3               !reference volume
+      EF_ref    = qe*n_ref*L_ref/eps_0   !reference electric field
+      Volt_ref  = EF_ref*L_ref           !reference voltage
 
     END SUBROUTINE readReference
 
@@ -94,6 +107,7 @@ MODULE moduleInput
       USE moduleRefParam
       USE moduleErrors
       USE moduleCaseParam
+      USE moduleSolver
       USE json_module
       IMPLICIT NONE
 
@@ -101,6 +115,7 @@ MODULE moduleInput
       LOGICAL:: found
       CHARACTER(:), ALLOCATABLE:: object
       REAL(8):: time !simulation time in [t]
+      CHARACTER(:), ALLOCATABLE:: solverType
 
       object = 'case'
       !Time parameters
@@ -110,6 +125,21 @@ MODULE moduleInput
       CALL config%get(object // '.time', time, found)
       IF (.NOT. found) CALL criticalError('Required parameter time not found','readCase')
 
+      CALL config%get(object // '.solver', solverType, found)
+      IF (.NOT. found) CALL criticalError('Required parameter solver not found','readCase')
+
+      !Allocate the type of solver
+      SELECT CASE(solverType)
+      CASE('2DCylNeutral')
+        !Allocate the solver for neutral pusher 2D Cylindrical
+        ALLOCATE(solver, source=solverCylNeutral())
+
+      CASE('2DCylCharged')
+        !Allocate the solver for charged pusher 2D Cylindrical
+        ALLOCATE(solver, source=solverCylCharged())
+
+      END SELECT
+
       !Convert simulation time to number of iterations
       tmax = INT(time/tau)
 
@@ -149,6 +179,7 @@ MODULE moduleInput
 
       CALL config%get(object // '.cpuTime', timeOutput, found)
       CALL config%get(object // '.numColl', collOutput, found)
+      CALL config%get(object // '.EMField', emOutput,   found)
 
     END SUBROUTINE readOutput
 
@@ -164,12 +195,12 @@ MODULE moduleInput
       CHARACTER(2):: iString
       CHARACTER(:), ALLOCATABLE:: object
       CHARACTER(:), ALLOCATABLE:: speciesType
-      REAL(8):: mass
+      REAL(8):: mass, charge
       LOGICAL:: found
       INTEGER:: i
 
       !Gets the number of species
-      CALL config%info('species', n_children = nSpecies)
+      CALL config%info('species', found, n_children = nSpecies)
       !Zero species means critical error
       IF (nSpecies == 0) CALL criticalError("No species found", "configRead")
 
@@ -180,21 +211,28 @@ MODULE moduleInput
         WRITE(iString, '(I2)') i
         object = 'species(' // TRIM(iString) // ')'
         CALL config%get(object // '.type', speciesType, found)
+        CALL config%get(object // '.mass', mass,        found)
+        mass = mass/m_ref
 
+        !Allocate species depending on type and assign specific parameters
         SELECT CASE(speciesType)
         CASE ("neutral")
           ALLOCATE(species(i)%obj, source=speciesNeutral())
-          !TODO: move to subroutine
-          CALL config%get(object // '.name',   species(i)%obj%name,        found)
-          CALL config%get(object // '.mass',   mass,                       found)
-          CALL config%get(object // '.weight', species(i)%obj%weight,      found)
-          species(i)%obj%pt = i
-          species(i)%obj%m  = mass/m_ref
+
+        CASE ("charged")
+          CALL config%get(object // '.charge', charge,      found)
+          ALLOCATE(species(i)%obj, source=speciesCharged(q  = charge, &
+                                                         qm = charge/mass))
 
         CASE DEFAULT
           CALL warningError("Species " // speciesType // " not supported yet")
 
         END SELECT
+        !Assign shared parameters for all species
+        CALL config%get(object // '.name',   species(i)%obj%name,   found)
+        CALL config%get(object // '.weight', species(i)%obj%weight, found)
+        species(i)%obj%pt = i
+        species(i)%obj%m  = mass
 
       END DO
 
@@ -217,29 +255,30 @@ MODULE moduleInput
       CHARACTER(:), ALLOCATABLE:: species_i, species_j
       CHARACTER(:), ALLOCATABLE:: crossSecFile
       CHARACTER(:), ALLOCATABLE:: crossSecFilePath
+      LOGICAL:: found
       INTEGER:: nInteractions, nCollisions
       INTEGER:: i, k, ij
       INTEGER:: pt_i, pt_j
 
       CALL initInteractionMatrix(interactionMatrix)
 
-      CALL config%get('interactions.folderCollisions', pathCollisions)
+      CALL config%get('interactions.folderCollisions', pathCollisions, found)
 
-      CALL config%info('interactions.collisions', n_children = nInteractions)
+      CALL config%info('interactions.collisions', found, n_children = nInteractions)
       DO i = 1, nInteractions
         WRITE(iString, '(I2)') i
         object = 'interactions.collisions(' // TRIM(iString) // ')'
-        CALL config%get(object // '.species_i', species_i)
+        CALL config%get(object // '.species_i', species_i, found)
         pt_i = speciesName2Index(species_i)
-        CALL config%get(object // '.species_j', species_j)
+        CALL config%get(object // '.species_j', species_j, found)
         pt_j = speciesName2Index(species_j)
-        CALL config%info(object // '.crossSections', n_children = nCollisions)
+        CALL config%info(object // '.crossSections', found, n_children = nCollisions)
         ij = interactionIndex(pt_i,pt_j)
         CALL interactionMatrix(ij)%init(nCollisions)
 
         DO k = 1, nCollisions
           WRITE (kString, '(I2)') k
-          CALL config%get(object // '.crossSections(' // TRIM(kString)// ')', crossSecFile)
+          CALL config%get(object // '.crossSections(' // TRIM(kString)// ')', crossSecFile, found)
           crossSecFilePath = pathCollisions // crossSecFile
           CALL interactionMatrix(ij)%collisions(k)%obj%init(crossSecFilePath, species(pt_i)%obj%m, species(pt_j)%obj%m)
 
@@ -247,7 +286,6 @@ MODULE moduleInput
 
       END DO
 
-
     END SUBROUTINE readInteractions
 
     !Reads boundary conditions for the mesh
@@ -263,7 +301,7 @@ MODULE moduleInput
       LOGICAL:: found
       INTEGER:: i
 
-      CALL config%info('boundary', n_children = nBoundary)
+      CALL config%info('boundary', found, n_children = nBoundary)
       ALLOCATE(boundary(1:nBoundary))
       DO i = 1, nBoundary
         WRITE(istring, '(i2)') i
@@ -322,6 +360,85 @@ MODULE moduleInput
 
     END SUBROUTINE readGeometry
 
+    SUBROUTINE readEMBoundary(config)
+      USE moduleMesh
+      USE moduleOutput
+      USE moduleErrors
+      USE moduleEM
+      USE moduleRefParam
+      USE moduleSpecies
+      USE json_module
+      IMPLICIT NONE
+
+      TYPE(json_file), INTENT(inout):: config
+      CHARACTER(:), ALLOCATABLE:: object
+      LOGICAL:: found
+      CHARACTER(2):: istring
+      INTEGER:: i, e, s
+      INTEGER:: info
+      EXTERNAL:: dgetrf
+
+      CALL config%info('boundaryEM', found, n_children = nBoundaryEM)
+      
+      IF (found) ALLOCATE(boundEM(1:nBoundaryEM))
+      
+      DO i = 1, nBoundaryEM
+        WRITE(istring, '(i2)') i
+        object = 'boundaryEM(' // trim(istring) // ')'
+
+        CALL config%get(object // '.type', boundEM(i)%typeEM, found)
+
+        SELECT CASE(boundEM(i)%typeEM)
+        CASE ("dirichlet")
+          CALL config%get(object // '.potential', boundEM(i)%potential, found)
+          IF (.NOT. found) &
+            CALL criticalError('Required parameter "potential" for Dirichlet boundary condition not found', 'readEMBoundary')
+          boundEM(i)%potential = boundEM(i)%potential/Volt_ref
+
+          CALL config%get(object // '.physicalSurface', boundEM(i)%physicalSurface, found)
+          IF (.NOT. found) &
+            CALL criticalError('Required parameter "physicalSurface" for Dirichlet boundary condition not found', 'readEMBoundary')
+
+        CASE DEFAULT
+          CALL criticalError('Boundary type ' // boundEM(i)%typeEM // ' not yet supported', 'readEMBoundary')
+
+        END SELECT
+
+      END DO
+
+      ALLOCATE(qSpecies(1:nSpecies))
+      DO s = 1, nSpecies
+        SELECT TYPE(sp => species(s)%obj)
+        TYPE IS (speciesCharged)
+          qSpecies(s) = sp%q
+
+        CLASS DEFAULT
+          qSpecies(s) = 0.D0
+
+        END SELECT
+
+      END DO
+
+      !TODO: Improve this
+      IF (ALLOCATED(boundEM)) THEN
+        DO e = 1, mesh%numEdges
+          IF (ANY(mesh%edges(e)%obj%physicalSurface == boundEM(:)%physicalSurface)) THEN
+            DO i = 1, nBoundaryEM
+              IF (mesh%edges(e)%obj%physicalSurface == boundEM(i)%physicalSurface) THEN
+                CALL boundEM(i)%apply(mesh%edges(e)%obj)
+
+              END IF
+            END DO
+          END IF
+        END DO
+      END IF
+
+      !Compute the PLU factorization of K once boundary conditions have been read
+      CALL dgetrf(mesh%numNodes, mesh%numNodes, mesh%K, mesh%numNodes, mesh%IPIV, info)
+      IF (info /= 0) CALL criticalError('Factorization of K matrix failed', 'readEMBoundary')
+
+    END SUBROUTINE readEMBoundary
+
     !Reads the injection of particles from the boundaries
     SUBROUTINE readInject(config)
       USE moduleSpecies
@@ -340,10 +457,11 @@ MODULE moduleInput
       REAL(8):: v
       REAL(8), ALLOCATABLE:: T(:), normal(:)
       REAL(8):: flow
+      CHARACTER(:), ALLOCATABLE:: units
       INTEGER:: physicalSurface
       INTEGER:: pt
 
-      CALL config%info('inject', n_children = nInject)
+      CALL config%info('inject', found, n_children = nInject)
       ALLOCATE(inject(1:nInject))
       nPartInj = 0
       DO i = 1, nInject
@@ -359,11 +477,13 @@ MODULE moduleInput
         CALL config%get(object // '.T',               T,               found)
         CALL config%get(object // '.n',               normal,          found)
         CALL config%get(object // '.flow',            flow,            found)
+        CALL config%get(object // '.units',           units,           found)
         CALL config%get(object // '.physicalSurface', physicalSurface, found)
 
-        CALL inject(i)%init(i, v, normal, T, flow, pt, physicalSurface)
+        CALL inject(i)%init(i, v, normal, T, flow, units, pt, physicalSurface)
 
       END DO
+      PRINT *, inject(:)%nParticles
 
       !Allocate array for injected particles
       IF (nPartInj > 0) THEN
diff --git a/src/modules/moduleMesh.f95 b/src/modules/moduleMesh.f95
index 85048ba..0b6e37f 100644
--- a/src/modules/moduleMesh.f95
+++ b/src/modules/moduleMesh.f95
@@ -2,7 +2,6 @@
 MODULE moduleMesh
   USE moduleList
   USE moduleOutput
-  USE OMP_LIB
   IMPLICIT NONE
 
   !Parent of Node element
@@ -13,6 +12,7 @@ MODULE moduleMesh
     REAL(8):: v = 0.D0
     !Output values
     TYPE(outputNode), ALLOCATABLE:: output(:)
+    TYPE(emNode):: emData
     CONTAINS
       PROCEDURE(initNode_interface), DEFERRED, PASS:: init
       PROCEDURE(getCoord_interface), DEFERRED, PASS:: getCoordinates
@@ -58,15 +58,17 @@ MODULE moduleMesh
     !id for boundary condition
     INTEGER:: bt = 0
     CONTAINS
-      PROCEDURE(initEdge_interface), DEFERRED, PASS:: init
-      PROCEDURE(boundary_interface), DEFERRED, PASS:: fBoundary
-      PROCEDURE(randPos_interface), DEFERRED, PASS::  randPos
+      PROCEDURE(initEdge_interface),     DEFERRED, PASS:: init
+      PROCEDURE(boundary_interface),     DEFERRED, PASS:: fBoundary
+      PROCEDURE(getNodesEdge_interface), DEFERRED, PASS:: getNodes
+      PROCEDURE(randPos_interface),      DEFERRED, PASS:: randPos
 
   END TYPE meshEdge
 
   ABSTRACT INTERFACE
     SUBROUTINE initEdge_interface(self, n, p, bt, physicalSurface)
       IMPORT:: meshEdge
+
       CLASS(meshEdge), INTENT(out):: self
       INTEGER, INTENT(in):: n
       INTEGER, INTENT(in):: p(:)
@@ -84,10 +86,16 @@ MODULE moduleMesh
 
     END SUBROUTINE
 
-    PURE FUNCTION randPos_interface(self, rnd) RESULT(r)
+    PURE FUNCTION getNodesEdge_interface(self) RESULT(n)
+      IMPORT:: meshEdge
+      CLASS(meshEdge), INTENT(in):: self
+      INTEGER, ALLOCATABLE:: n(:)
+
+    END FUNCTION
+
+    FUNCTION randPos_interface(self) RESULT(r)
       IMPORT:: meshEdge
       CLASS(meshEdge), INTENT(in):: self
-      REAL(8), INTENT(in):: rnd
       REAL(8):: r(1:3)
 
     END FUNCTION randPos_interface
@@ -119,6 +127,9 @@ MODULE moduleMesh
     CONTAINS
       PROCEDURE(initVol_interface),     DEFERRED, PASS:: init
       PROCEDURE(scatter_interface),     DEFERRED, PASS:: scatter
+      PROCEDURE(gatherEF_interface),    DEFERRED, PASS:: gatherEF
+      PROCEDURE(getNodesVol_interface), DEFERRED, PASS:: getNodes
+      PROCEDURE(elemF_interface),       DEFERRED, PASS:: elemF
       PROCEDURE(collision_interface),   DEFERRED, PASS:: collision
       PROCEDURE(findCell_interface),    DEFERRED, PASS:: findCell
       PROCEDURE(resetOutput_interface), DEFERRED, PASS:: resetOutput
@@ -143,6 +154,32 @@ MODULE moduleMesh
 
     END SUBROUTINE scatter_interface
 
+    PURE FUNCTION gatherEF_interface(self, xi) RESULT(EF)
+      
+      IMPORT:: meshVol
+      CLASS(meshVol), INTENT(in):: self
+      REAL(8), INTENT(in):: xi(1:3)
+      REAL(8):: EF(1:3)
+
+    END FUNCTION gatherEF_interface
+
+    PURE FUNCTION getNodesVol_interface(self) RESULT(n)
+      
+      IMPORT:: meshVol
+      CLASS(meshVol), INTENT(in):: self
+      INTEGER, ALLOCATABLE:: n(:)
+
+    END FUNCTION getNodesVol_interface
+
+    PURE FUNCTION elemF_interface(self, source) RESULT(localF)
+      
+      IMPORT:: meshVol
+      CLASS(meshVol), INTENT(in):: self
+      REAL(8), INTENT(in):: source(1:)
+      REAL(8), ALLOCATABLE:: localF(:)
+
+    END FUNCTION elemF_interface
+
     SUBROUTINE collision_interface(self)
       IMPORT:: meshVol
       CLASS(meshVol), INTENT(inout):: self
@@ -184,13 +221,14 @@ MODULE moduleMesh
     TYPE(meshVolCont),  ALLOCATABLE:: vols(:)
     !Global stiffness matrix
     REAL(8), ALLOCATABLE, DIMENSION(:,:):: K 
-    !Global load vector
-    REAL(8), ALLOCATABLE, DIMENSION(:):: F 
+    !Permutation matrix for P L U factorization
+    INTEGER, ALLOCATABLE, DIMENSION(:,:):: IPIV
 
     CONTAINS
       PROCEDURE(readMesh_interface),    PASS, DEFERRED:: readMesh
       PROCEDURE(printOutput_interface), PASS, DEFERRED:: printOutput
       PROCEDURE(printColl_interface),   PASS, DEFERRED:: printColl
+      PROCEDURE(printEM_interface),     PASS, DEFERRED:: printEM
 
   END TYPE meshGeneric
 
@@ -213,7 +251,7 @@ MODULE moduleMesh
 
     END SUBROUTINE printOutput_interface
 
-    !Prints unmber of collisions
+    !Prints number of collisions
     SUBROUTINE printColl_interface(self, t)
       IMPORT meshGeneric
 
@@ -222,6 +260,14 @@ MODULE moduleMesh
 
     END SUBROUTINE printColl_interface
 
+    !Prints EM info
+    SUBROUTINE printEM_interface(self, t)
+      IMPORT meshGeneric
+
+      CLASS(meshGeneric), INTENT(in):: self
+      INTEGER, INTENT(in):: t
+
+    END SUBROUTINE printEM_interface
   END INTERFACE
 
   !Generic mesh
diff --git a/src/modules/moduleMeshCyl.f95 b/src/modules/moduleMeshCyl.f95
index 69c9bc6..ac912c8 100644
--- a/src/modules/moduleMeshCyl.f95
+++ b/src/modules/moduleMeshCyl.f95
@@ -25,8 +25,9 @@ MODULE moduleMeshCyl
     !Connectivity to nodes
     CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL()
     CONTAINS
-      PROCEDURE, PASS:: init    => initEdgeCyl
-      PROCEDURE, PASS:: randPos => randPosCyl
+      PROCEDURE, PASS:: init     => initEdgeCyl
+      PROCEDURE, PASS:: getNodes => getNodesCyl
+      PROCEDURE, PASS:: randPos  => randPosCyl
 
   END TYPE meshEdgeCyl
 
@@ -43,22 +44,20 @@ MODULE moduleMeshCyl
     CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL(), n3 => NULL(), n4 => NULL()
     !Connectivity to adjacent elements
     CLASS(*), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL(), e4 => NULL()
-    !Maximum collision rate
-    !TODO: Move to other more appropiate section
-    REAL(8):: phi(1:4) = 0.D0
-    REAL(8):: Ke(1:4,1:4) = 0.D0
     REAL(8):: arNodes(1:4) = 0.D0
 
     CONTAINS
       PROCEDURE, PASS::   init        => initVolQuadCyl
-      PROCEDURE, PASS::   locKe       => localKeQuad
+      PROCEDURE, PASS::   elemK       => elemKQuad
+      PROCEDURE, PASS::   elemF       => elemFQuad
       PROCEDURE, PASS::   detJac      => detJQuad
       PROCEDURE, PASS::   invJ        => invJQuad
       PROCEDURE, PASS::   area        => areaQuad
       PROCEDURE, NOPASS:: weight      => weightQuad
       PROCEDURE, NOPASS:: inside      => insideQuad
-      PROCEDURE, PASS::   dVal        => dValQuad
       PROCEDURE, PASS::   scatter     => scatterQuad
+      PROCEDURE, PASS::   gatherEF    => gatherEFQuad
+      PROCEDURE, PASS::   getNodes    => getNodesQuad
       PROCEDURE, PASS::   phy2log     => phy2logQuad
       PROCEDURE, PASS::   findCell    => findCellCylQuad
       PROCEDURE, PASS::   resetOutput => resetOutputQuad
@@ -78,8 +77,8 @@ MODULE moduleMeshCyl
       REAL(8), INTENT(in):: r(1:3)
 
       self%n = n
-      self%r = r(1)/L_ref
       self%z = r(2)/L_ref
+      self%r = r(1)/L_ref
       !Node volume, to be determined in mesh
       self%v = 0.D0
 
@@ -128,13 +127,25 @@ MODULE moduleMeshCyl
 
     END SUBROUTINE initEdgeCyl
 
+    PURE FUNCTION getNodesCyl(self) RESULT(n)
+      IMPLICIT NONE
 
-    PURE FUNCTION randPosCyl(self, rnd) RESULT(r)
       CLASS(meshEdgeCyl), INTENT(in):: self
-      REAL(8), INTENT(in):: rnd
+      INTEGER, ALLOCATABLE:: n(:)
+
+      ALLOCATE(n(1:2))
+      n = (/self%n1%n, self%n2%n /)
+
+    END FUNCTION getNodesCyl
+
+    FUNCTION randPosCyl(self) RESULT(r)
+      CLASS(meshEdgeCyl), INTENT(in):: self
+      REAL(8):: rnd
       REAL(8):: r(1:3)
       REAL(8):: p1(1:2), p2(1:2)
 
+      rnd = RAND()
+
       p1 = (/self%z(1), self%r(1) /)
       p2 = (/self%z(2), self%r(2) /)
       r(1:2) = (1.D0 - rnd)*p1 + rnd*p2
@@ -173,15 +184,14 @@ MODULE moduleMeshCyl
       self%n3%v = self%n3%v + self%arNodes(3)
       self%n4%v = self%n4%v + self%arNodes(4)
 
-      CALL self%locKe()
-
       self%sigmaVrelMax = sigma_ref/L_ref**2
 
       CALL OMP_INIT_LOCK(self%lock)
 
     END SUBROUTINE initVolQuadCyl
 
-    FUNCTION fpsi(xi1,xi2) RESULT(psi) 
+    !TODO: MAKE ELEMENT DEPENDENT
+    PURE FUNCTION fpsi(xi1,xi2) RESULT(psi) 
       IMPLICIT NONE
 
       REAL(8),INTENT(in):: xi1,xi2
@@ -195,7 +205,8 @@ MODULE moduleMeshCyl
     END FUNCTION fpsi
 
     !Derivative element function (xi1)
-    FUNCTION dpsiXi1(xi2) RESULT(psiXi1) 
+    !TODO: MAKE ELEMENT DEPENDENT
+    PURE FUNCTION dpsiXi1(xi2) RESULT(psiXi1) 
       IMPLICIT NONE
 
       REAL(8),INTENT(in):: xi2
@@ -209,7 +220,8 @@ MODULE moduleMeshCyl
     END FUNCTION dpsiXi1
 
     !Derivative element function (xi2)
-    FUNCTION dpsiXi2(xi1) RESULT(psiXi2) 
+    !TODO: MAKE ELEMENT DEPENDENT
+    PURE FUNCTION dpsiXi2(xi1) RESULT(psiXi2) 
       IMPLICIT NONE
 
       REAL(8),INTENT(in):: xi1
@@ -223,12 +235,13 @@ MODULE moduleMeshCyl
     END FUNCTION dpsiXi2
 
     !Transforms physical coordinates to element coordinates
-    FUNCTION phy2logQuad(this,r) RESULT(xN) 
+    PURE FUNCTION phy2logQuad(self,r) RESULT(xN) 
       IMPLICIT NONE
 
-      CLASS(meshVolCylQuad):: this
-      REAL(8):: r(1:2)
-      REAL(8):: xN(1:2), xO(1:2), dPsi(1:2,1:4), detJ, j(1:2,1:2), psi(1:4), f(1:2)
+      CLASS(meshVolCylQuad), INTENT(in):: self
+      REAL(8), INTENT(in):: r(1:3)
+      REAL(8):: xN(1:3)
+      REAL(8):: xO(1:3), dPsi(1:2,1:4), detJ, j(1:2,1:2), psi(1:4), f(1:2)
       REAL(8):: conv
 
       !Iterative newton method to transform coordinates
@@ -238,12 +251,12 @@ MODULE moduleMeshCyl
       DO WHILE(conv>1.D-4)
         dPsi(1,:) = dpsiXi1(xO(2))
         dPsi(2,:) = dpsiXi2(xO(1))
-        j    = this%invJ(xO(1),xO(2),dPsi)
-        psi = fpsi(xO(1), xO(2))
-        f(1) = DOT_PRODUCT(psi,this%z)-r(1)
-        f(2) = DOT_PRODUCT(psi,this%r)-r(2)
-        detJ = this%detJac(xO(1),xO(2),dPsi)
-        xN=xO - MATMUL(j, f)/detJ
+        j    = self%invJ(xO, dPsi)
+        psi  = fpsi(xO(1), xO(2))
+        f(1) = DOT_PRODUCT(psi,self%z)-r(1)
+        f(2) = DOT_PRODUCT(psi,self%r)-r(2)
+        detJ = self%detJac(xO,dPsi)
+        xN(1:2)=xO(1:2) - MATMUL(j, f)/detJ
         conv=MAXVAL(DABS(xN-xO),1)
         xO=xN
 
@@ -252,106 +265,139 @@ MODULE moduleMeshCyl
     END FUNCTION phy2logQuad
 
     !Computes element local stiffness matrix
-    SUBROUTINE localKeQuad(this) 
+    PURE FUNCTION elemKQuad(self) RESULT(ke)
       USE moduleConstParam
       IMPLICIT NONE
 
-      CLASS(meshVolCylQuad):: this
-      REAL(8):: r, xi1, xi2, dPsi(1:2,1:4)
+      CLASS(meshVolCylQuad), INTENT(in):: self
+      REAL(8):: r, xi(1:3), dPsi(1:2,1:4)
+      REAL(8):: ke(1:4,1:4)
       REAL(8):: j(1:2,1:2), detJ
       INTEGER:: l, m
 
-      this%Ke=0.D0
+      ke=0.D0
       !Start 2D Gauss Quad Integral
       DO l=1, 3
+        xi(1)     = cor(l)
+        dPsi(2,:) = dpsiXi2(xi(1))
         DO m = 1, 3
-          xi1 = cor(l)
-          xi2 = cor(m)
-          dPsi(1,:) = dpsiXi1(xi2)
-          dPsi(2,:) = dpsiXi2(xi1)
-          detJ = this%detJac(xi1,xi2,dPsi)
-          j    = this%invJ(xi1,xi2,dPsi)
-          r    = DOT_PRODUCT(fpsi(xi1,xi2),this%r)
-          this%Ke = this%Ke + MATMUL(TRANSPOSE(MATMUL(j,dPsi)),MATMUL(j,dPsi))*r*w(l)*w(m)/detJ
+          xi(2)     = cor(m)
+          dPsi(1,:) = dpsiXi1(xi(2))
+          detJ      = self%detJac(xi,dPsi)
+          j         = self%invJ(xi,dPsi)
+          r         = DOT_PRODUCT(fpsi(xi(1),xi(2)),self%r)
+          ke        = ke + MATMUL(TRANSPOSE(MATMUL(j,dPsi)),MATMUL(j,dPsi))*r*w(l)*w(m)/detJ
 
         END DO
       END DO
-      this%Ke = this%Ke*2.D0*PI
+      ke = ke*2.D0*PI
 
-    END SUBROUTINE localKeQuad
+    END FUNCTION elemKQuad
 
-    !Computes element Jacobian determinant
-    FUNCTION detJQuad(this,xi1,xi2,dPsi_in) RESULT(dJ) 
+    !Computes the local source vector for a force f
+    PURE FUNCTION elemFQuad(self, source) RESULT(localF)
+      USE moduleConstParam
       IMPLICIT NONE
 
-      REAL(8),OPTIONAL:: dPsi_in(1:2,1:4)
+      CLASS(meshVolCylQuad), INTENT(in):: self
+      REAL(8), INTENT(in):: source(1:)
+      REAL(8), ALLOCATABLE:: localF(:)
+      REAL(8):: r, xi(1:3), psi(1:4)
+      REAL(8):: detJ, f
+      INTEGER:: l, m
+
+      ALLOCATE(localF(1:4))
+      localF = 0.D0
+      xi = 0.D0
+      DO l=1, 3
+        DO m = 1, 3
+          xi(1) = cor(l)
+          xi(2) = cor(m)
+          detJ = self%detJac(xi)
+          psi = fpsi(xi(1),xi(2))
+          r = DOT_PRODUCT(psi,self%r)
+          f = DOT_PRODUCT(psi,source)
+          localF = localF + r*f*psi*w(l)*w(m)*detJ
+
+         END DO
+      END DO
+      localF = localF*2.D0*PI
+
+    END FUNCTION elemFQuad
+
+    !Computes element Jacobian determinant
+    PURE FUNCTION detJQuad(self,xi,dPsi_in) RESULT(dJ) 
+      IMPLICIT NONE
+
+      CLASS(meshVolCylQuad), INTENT(in):: self
+      REAL(8), INTENT(in):: xi(1:3)
+      REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:2,1:4)
+      REAL(8):: dJ
       REAL(8):: dPsi(1:2,1:4)
       REAL(8):: dz(1:2), dr(1:2)
-      REAL(8):: xi1,xi2, dJ
-      CLASS(meshVolCylQuad):: this
 
       IF(PRESENT(dPsi_in)) THEN
         dPsi=dPsi_in
 
       ELSE
-        dPsi(1,:) = dpsiXi1(xi2)
-        dPsi(2,:) = dpsiXi2(xi1)
+        dPsi(1,:) = dpsiXi1(xi(2))
+        dPsi(2,:) = dpsiXi2(xi(1))
 
       END IF
 
-      dz(1) = DOT_PRODUCT(dPsi(1,:),this%z)
-      dz(2) = DOT_PRODUCT(dPsi(2,:),this%z)
-      dr(1) = DOT_PRODUCT(dPsi(1,:),this%r)
-      dr(2) = DOT_PRODUCT(dPsi(2,:),this%r)
+      dz(1) = DOT_PRODUCT(dPsi(1,:),self%z)
+      dz(2) = DOT_PRODUCT(dPsi(2,:),self%z)
+      dr(1) = DOT_PRODUCT(dPsi(1,:),self%r)
+      dr(2) = DOT_PRODUCT(dPsi(2,:),self%r)
       dJ = dz(1)*dr(2)-dz(2)*dr(1)
 
     END FUNCTION detJQuad
 
-    FUNCTION invJQuad(this,xi1,xi2,dPsi_in) RESULT(j) !Computes element Jacobian inverse matrix (without determinant)
+    PURE FUNCTION invJQuad(self,xi,dPsi_in) RESULT(j) !Computes element Jacobian inverse matrix (without determinant)
       IMPLICIT NONE
 
-      REAL(8),OPTIONAL:: dpsi_in(1:2,1:4)
+      CLASS(meshVolCylQuad), INTENT(in):: self
+      REAL(8), INTENT(in):: xi(1:3)
+      REAL(8), INTENT(in), OPTIONAL:: dpsi_in(1:2,1:4)
       REAL(8):: dPsi(1:2,1:4)
       REAL(8):: dz(1:2), dr(1:2)
-      REAL(8):: xi1,xi2, j(1:2,1:2)
-      CLASS(meshVolCylQuad):: this
+      REAL(8):: j(1:2,1:2)
 
       IF(PRESENT(dPsi_in)) THEN
         dPsi=dPsi_in
 
       ELSE
-        dPsi(1,:) = dpsiXi1(xi2)
-        dPsi(2,:) = dpsiXi2(xi1)
+        dPsi(1,:) = dpsiXi1(xi(2))
+        dPsi(2,:) = dpsiXi2(xi(1))
 
       END IF
 
-      dz(1) = DOT_PRODUCT(dPsi(1,:),this%z)
-      dz(2) = DOT_PRODUCT(dPsi(2,:),this%z)
-      dr(1) = DOT_PRODUCT(dPsi(1,:),this%r)
-      dr(2) = DOT_PRODUCT(dPsi(2,:),this%r)
+      dz(1) = DOT_PRODUCT(dPsi(1,:),self%z)
+      dz(2) = DOT_PRODUCT(dPsi(2,:),self%z)
+      dr(1) = DOT_PRODUCT(dPsi(1,:),self%r)
+      dr(2) = DOT_PRODUCT(dPsi(2,:),self%r)
       j(1,:) = (/  dr(2), -dz(2) /)
       j(2,:) = (/ -dr(1),  dz(1) /)
 
     END FUNCTION invJQuad
 
-    SUBROUTINE areaQuad(this)!Computes element area
+    SUBROUTINE areaQuad(self)!Computes element area
       USE moduleConstParam
       IMPLICIT NONE
 
-      CLASS(meshVolCylQuad):: this
-      REAL(8):: r, xi1, xi2
+      CLASS(meshVolCylQuad):: self
+      REAL(8):: r, xi(1:3)
       REAL(8):: detJ, fpsi0(1:4)
 
-      this%volume  = 0.D0
-      this%arNodes = 0.D0
+      self%volume  = 0.D0
+      self%arNodes = 0.D0
       !2D 1 point Gauss Quad Integral
-      xi1 = 0.D0
-      xi2 = 0.D0
-      detJ = this%detJac(xi1,xi2)*8.D0*PI !4*2*pi
-      fpsi0 = fpsi(xi1,xi2)
-      r = DOT_PRODUCT(fpsi0,this%r)
-      this%volume  =       r*detJ
-      this%arNodes = fpsi0*r*detJ
+      xi = 0.D0
+      detJ = self%detJac(xi)*8.D0*PI !4*2*pi
+      fpsi0 = fpsi(xi(1),xi(2))
+      r = DOT_PRODUCT(fpsi0,self%r)
+      self%volume  =       r*detJ
+      self%arNodes = fpsi0*r*detJ
 
     END SUBROUTINE areaQuad
 
@@ -376,24 +422,6 @@ MODULE moduleMeshCyl
 
     END FUNCTION insideQuad
 
-    FUNCTION dValQuad(this,xi1,xi2) RESULT(EF)!Computes electric field in xi1,xi2
-      IMPLICIT NONE
-
-      CLASS(meshVolCylQuad):: this
-      REAL(8):: xi1, xi2, dPsi(1:2,1:4)
-      REAL(8):: j(1:2,1:2), detJ
-      REAL(8):: EF(1:3)
-
-      dPsi(1,:) = dpsiXi1(xi2)
-      dPsi(2,:) = dpsiXi2(xi1)
-      detJ = this%detJac(xi1,xi2,dPsi)
-      j    = this%invJ(xi1,xi2,dPsi)
-      EF(1) = -DOT_PRODUCT(dPsi(1,:), this%phi)*j(2,2)/detJ
-      EF(2) = -DOT_PRODUCT(dPsi(2,:), this%phi)*j(1,1)/detJ
-      EF(3) =  0.D0
-
-    END FUNCTION dValQuad
-
     SUBROUTINE scatterQuad(self, part)
       USE moduleOutput
       USE moduleSpecies
@@ -428,29 +456,65 @@ MODULE moduleMeshCyl
       vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(4)*part%v(:)
       vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(4)*tensorS
 
-
     END SUBROUTINE scatterQuad
 
+    PURE FUNCTION gatherEFQuad(self,xi) RESULT(EF)!Computes electric field inside element
+      IMPLICIT NONE
+
+      CLASS(meshVolCylQuad), INTENT(in):: self
+      REAL(8), INTENT(in):: xi(1:3)
+      REAL(8):: dPsi(1:2,1:4)
+      REAL(8):: j(1:2,1:2), detJ
+      REAL(8):: phi(1:4)
+      REAL(8):: EF(1:3)
+
+      phi = (/self%n1%emData%phi, &
+              self%n2%emData%phi, &
+              self%n3%emData%phi, &
+              self%n4%emData%phi /)
+
+      dPsi(1,:) = dpsiXi1(xi(2))
+      dPsi(2,:) = dpsiXi2(xi(1))
+      detJ = self%detJac(xi,dPsi)
+      j    = self%invJ(xi,dPsi)
+      EF(1) = -DOT_PRODUCT(dPsi(1,:), phi)*j(2,2)/detJ
+      EF(2) = -DOT_PRODUCT(dPsi(2,:), phi)*j(1,1)/detJ
+      EF(3) =  0.D0
+
+    END FUNCTION gatherEFQuad
+
+    PURE FUNCTION getNodesQuad(self) RESULT(n)
+      IMPLICIT NONE
+
+      CLASS(meshVolCylQuad), INTENT(in):: self
+      INTEGER, ALLOCATABLE:: n(:)
+
+      ALLOCATE(n(1:4))
+      n = (/self%n1%n, self%n2%n, self%n3%n, self%n4%n /)
+
+    END FUNCTION getNodesQuad
+
     RECURSIVE SUBROUTINE findCellCylQuad(self, part, oldCell)
       USE moduleSpecies
+      USE OMP_LIB
       IMPLICIT NONE
 
       CLASS(meshVolCylQuad), INTENT(inout):: self
       CLASS(meshVol), OPTIONAL, INTENT(in):: oldCell
       CLASS(particle), INTENT(inout), TARGET:: part
-      REAL(8):: xLog(1:2)
+      REAL(8):: xLog(1:3)
       REAL(8):: xLogArray(1:4)
       CLASS(*), POINTER:: nextElement
       INTEGER:: nextInt
 
-      xLog = self%phy2log(part%r(1:2))
+      xLog = self%phy2log(part%r)
+      !Checks if particle is inside of current cell
       IF (self%inside(xLog(1), xLog(2))) THEN
-        !Checks if particle is inside of current cell
-        IF (PRESENT(oldCell)) THEN
-          !If oldCell, recalculate particle weight, as particle has entered a new cell.
-          part%weight = part%weight*oldCell%volume/self%volume
-
-        END IF
+        ! IF (PRESENT(oldCell)) THEN
+        !   !If oldCell, recalculate particle weight, as particle has entered a new cell.
+        !   part%weight = part%weight*oldCell%volume/self%volume
+        !
+        ! END IF
         part%e_p  = self%n
         part%xLog = xLog
         part%n_in = .TRUE.
diff --git a/src/modules/moduleMeshCylBoundary.f95 b/src/modules/moduleMeshCylBoundary.f95
index 4298193..7c29115 100644
--- a/src/modules/moduleMeshCylBoundary.f95
+++ b/src/modules/moduleMeshCylBoundary.f95
@@ -74,8 +74,6 @@ MODULE moduleMeshCylBoundary
       CLASS(meshEdgeCylAxis), INTENT(inout):: self
       CLASS(particle), INTENT(inout):: part
 
-      !Do to the 2D Cyl pusher, this function should never be called.
-
     END SUBROUTINE symmetryAxis
 
 
diff --git a/src/modules/moduleMeshCylRead.f95 b/src/modules/moduleMeshCylRead.f95
index df6e51b..9a0f5ff 100644
--- a/src/modules/moduleMeshCylRead.f95
+++ b/src/modules/moduleMeshCylRead.f95
@@ -8,6 +8,7 @@ MODULE moduleMeshCylRead
       PROCEDURE, PASS:: readMesh => readMeshCyl
       PROCEDURE, PASS:: printOutput => printOutputCyl
       PROCEDURE, PASS:: printColl   => printCollisionsCyl
+      PROCEDURE, PASS:: printEM     => printEMCyl
 
   END TYPE
 
@@ -41,8 +42,10 @@ MODULE moduleMeshCylRead
       READ(10, *) self%numNodes
       !Allocate required matrices and vectors
       ALLOCATE(self%nodes(1:self%numNodes))
-      ! ALLOCATE(self%K(1:numNodes,1:numNodes))
-      ! ALLOCATE(self%F(1:numNodes))
+      ALLOCATE(self%K(1:self%numNodes,1:self%numNodes))
+      ALLOCATE(self%IPIV(1:self%numNodes,1:self%numNodes))
+      self%K = 0.D0
+      self%IPIV = 0
       !Read nodes cartesian coordinates (x=z, y=r, z=null)
       DO e=1, self%numNodes
         READ(10, *) n, z, r
@@ -123,24 +126,10 @@ MODULE moduleMeshCylRead
 
         END DO
 
-      END DO
+        !Constructs the global K matrix
+        CALL constructGlobalK(self%K, self%vols(e)%obj)
 
-      ! !Compute global stiffness matrix
-      ! GlobalK=0.D0
-      ! DO e=1, numElem
-        ! DO i=1, 4
-          ! DO j=1, 4
-            ! GlobalK(elems(e)%p(i),elems(e)%p(j)) = GlobalK(elems(e)%p(i),elems(e)%p(j)) + elems(e)%Ke(i,j)
-          ! END DO
-        ! END DO
-      ! END DO
-      ! ! Apply Dirichlet boundary conditions to GlobalK
-      ! DO n=1, numNodes
-        ! IF (nodes(n)%bound == 1) THEN
-          ! GlobalK(n,:)=0.D0
-          ! GlobalK(n,n)=1.D0
-        ! END IF
-      ! END DO
+      END DO
 
     END SUBROUTINE
 
@@ -298,6 +287,41 @@ MODULE moduleMeshCylRead
 
     END SUBROUTINE connectedQuadEdge
 
+    SUBROUTINE constructGlobalK(K, elem)
+      IMPLICIT NONE
+
+      REAL(8), INTENT(inout):: K(1:,1:)
+      CLASS(meshVol), INTENT(in):: elem
+      REAL(8), ALLOCATABLE:: localK(:,:)
+      INTEGER:: nNodes, i, j
+      INTEGER, ALLOCATABLE:: n(:)
+
+      SELECT TYPE(elem)
+      TYPE IS(meshVolCylQuad)
+        nNodes = 4
+        ALLOCATE(localK(1:nNodes,1:nNodes))
+        localK = elem%elemK()
+        ALLOCATE(n(1:nNodes))
+        n = (/ elem%n1%n, elem%n2%n, &
+               elem%n3%n, elem%n4%n /)
+
+      CLASS DEFAULT
+        nNodes = 0
+        ALLOCATE(localK(1:1, 1:1))
+        localK = 0.D0
+        ALLOCATE(n(1:1))
+        n = 0
+
+      END SELECT
+      
+      DO i = 1, nNodes
+        DO j = 1, nNodes
+          K(n(i), n(j)) = K(n(i), n(j)) + localK(i, j)
+        END DO
+      END DO
+
+    END SUBROUTINE constructGlobalK
+
     SUBROUTINE printOutputCyl(self, t)
       USE moduleRefParam
       USE moduleSpecies
@@ -425,4 +449,60 @@ MODULE moduleMeshCylRead
 
     END SUBROUTINE printCollisionsCyl
 
+    SUBROUTINE printEMCyl(self, t)
+      USE moduleRefParam
+      USE moduleCaseParam
+      USE moduleOutput
+      IMPLICIT NONE
+
+      CLASS(meshCylGeneric), INTENT(in):: self
+      INTEGER, INTENT(in):: t
+      INTEGER:: n, e
+      REAL(8):: time
+      CHARACTER(:), ALLOCATABLE:: fileName
+      CHARACTER (LEN=6):: tstring !TODO: Review to allow any number of iterations
+
+      IF (emOutput) THEN
+        time = DBLE(t)*tau*ti_ref
+        WRITE(tstring, '(I6.6)') t
+
+        fileName='OUTPUT_' // tstring// '_EMField.msh'
+        WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
+        OPEN (20, file = path // folder // '/' //  fileName)
+        WRITE(20, "(A)") '$MeshFormat'
+        WRITE(20, "(A)") '2.2 0 8'
+        WRITE(20, "(A)") '$EndMeshFormat'
+        WRITE(20, "(A)") '$NodeData'
+        WRITE(20, "(A)") '1'
+        WRITE(20, "(A)") '"Potential (V)"'
+        WRITE(20, *) 1
+        WRITE(20, *) time
+        WRITE(20, *) 3
+        WRITE(20, *) t
+        WRITE(20, *) 1
+        WRITE(20, *) self%numNodes
+        DO n=1, self%numNodes
+          WRITE(20, *) n, self%nodes(n)%obj%emData%phi*Volt_ref
+        END DO
+        WRITE(20, "(A)") '$EndNodeData'
+        
+        WRITE(20, "(A)") '$ElementData'
+        WRITE(20, "(A)") '1'
+        WRITE(20, "(A)") '"Electric Field (V/m)"'
+        WRITE(20, *) 1
+        WRITE(20, *) time
+        WRITE(20, *) 3
+        WRITE(20, *) t
+        WRITE(20, *) 3
+        WRITE(20, *) self%numVols
+        DO e=1, self%numVols
+          WRITE(20, *) e+self%numEdges, self%vols(e)%obj%gatherEF((/0.D0, 0.D0, 0.D0/))*EF_ref
+        END DO
+        WRITE(20, "(A)") '$EndElementData'
+        CLOSE(20)
+
+      END IF
+
+    END SUBROUTINE printEMCyl
+  
 END MODULE moduleMeshCylRead
diff --git a/src/modules/moduleOutput.f95 b/src/modules/moduleOutput.f95
index dfe3808..fe1df6f 100644
--- a/src/modules/moduleOutput.f95
+++ b/src/modules/moduleOutput.f95
@@ -7,6 +7,13 @@ MODULE moduleOutput
 
   END TYPE
 
+  !Type for EM data in node
+  TYPE emNode
+    CHARACTER(:), ALLOCATABLE:: type
+    REAL(8):: phi
+
+  END TYPE emNode
+
   !Output in dimensional units to print
   TYPE outputFormat
     REAL(8):: density, velocity(1:3), pressure, temperature
@@ -18,6 +25,7 @@ MODULE moduleOutput
   INTEGER:: triggerOutput, counterOutput
   LOGICAL:: timeOutput = .FALSE.
   LOGICAL:: collOutput = .FALSE.
+  LOGICAL:: emOutput   = .FALSE.
 
   CONTAINS
     FUNCTION outerProduct(a,b) RESULT(s)
@@ -42,6 +50,7 @@ MODULE moduleOutput
     END FUNCTION tensorTrace
 
     SUBROUTINE calculateOutput(rawValues, formatValues, nodeVol, speciesIn)
+      USE moduleConstParam
       USE moduleRefParam
       USE moduleSpecies
       IMPLICIT NONE
diff --git a/src/modules/moduleRefParam.f95 b/src/modules/moduleRefParam.f95
index 152f103..9d2b8d9 100644
--- a/src/modules/moduleRefParam.f95
+++ b/src/modules/moduleRefParam.f95
@@ -1,10 +1,9 @@
 !Reference parameters
 MODULE moduleRefParam
-  USE moduleConstParam
   !Parameters that define the problem (inputs)
-  REAL(8):: n_ref, m_ref, T_ref, r_ref, sigma_ref!, q_ref=1.6022D-19
+  REAL(8):: n_ref, m_ref, T_ref, r_ref, debye_ref, sigma_ref
   !Reference parameters for non-dimensional problem
-  REAL(8):: L_ref, v_ref, ti_ref, Vol_ref
+  REAL(8):: L_ref, v_ref, ti_ref, Vol_ref, EF_ref, Volt_ref
 
 END MODULE moduleRefParam
 
diff --git a/src/modules/moduleSolver.f95 b/src/modules/moduleSolver.f95
index ec05cd7..21b3766 100644
--- a/src/modules/moduleSolver.f95
+++ b/src/modules/moduleSolver.f95
@@ -1,5 +1,44 @@
 MODULE moduleSolver
 
+  TYPE, PUBLIC, ABSTRACT:: solverGeneric
+    CONTAINS
+      PROCEDURE(push_interafece), DEFERRED, NOPASS:: push
+      PROCEDURE(solveEM_interface), DEFERRED, NOPASS:: solveEMField
+
+  END TYPE solverGeneric
+
+  ABSTRACT INTERFACE
+    !Push a particle
+    PURE SUBROUTINE push_interafece(part)
+      USE moduleSpecies
+
+      TYPE(particle), INTENT(inout):: part
+
+    END SUBROUTINE push_interafece
+
+    !Solve the electromagnetic field
+    SUBROUTINE solveEM_interface()
+    
+    END SUBROUTINE solveEM_interface
+
+  END INTERFACE
+
+  TYPE, PUBLIC, EXTENDS(solverGeneric):: solverCylNeutral
+    CONTAINS
+      PROCEDURE, NOPASS:: push => pushCylNeutral
+      PROCEDURE, NOPASS:: solveEMField => noEMField
+
+  END TYPE solverCylNeutral
+
+  TYPE, PUBLIC, EXTENDS(solverGeneric):: solverCylCharged
+    CONTAINS
+      PROCEDURE, NOPASS:: push => pushCylCharged
+      PROCEDURE, NOPASS:: solveEMField => solveElecField
+
+  END TYPE solverCylCharged
+
+  CLASS(solverGeneric), ALLOCATABLE:: solver
+
   CONTAINS
     SUBROUTINE doPushes()
       USE moduleSpecies
@@ -10,7 +49,9 @@ MODULE moduleSolver
 
       !$OMP DO
       DO n=1, nPartOld
-        CALL push(partOld(n))
+        !Push particle
+        CALL solver%push(partOld(n))
+        !Find cell in wich particle reside
         CALL mesh%vols(partOld(n)%e_p)%obj%findCell(partOld(n))
 
       END DO
@@ -18,16 +59,15 @@ MODULE moduleSolver
 
     END SUBROUTINE doPushes
 
-    !Push one particle. Boris pusher for 2D Cyl
-    !TODO: make general for any pusher
-    PURE SUBROUTINE push(part)
+    !Push one particle. Boris pusher for 2D Cyl Neutral particle
+    PURE SUBROUTINE pushCylNeutral(part)
       USE moduleSpecies
-      USE moduleMesh
       IMPLICIT NONE
+
       TYPE(particle), INTENT(inout):: part
       REAL(8):: v_p_oh_star(2:3)
       TYPE(particle):: part_temp
-      REAL(8):: x_new, y_new, r, sin_alpha, cos_alpha, alpha
+      REAL(8):: x_new, y_new, r, sin_alpha, cos_alpha
 
       part_temp = part
       !z
@@ -40,8 +80,6 @@ MODULE moduleSolver
       y_new =             v_p_oh_star(3)*tau
       r = DSQRT(x_new**2+y_new**2)
       part_temp%r(2) = r
-      alpha = 0.D0!ATAN2(y_new,x_new) !0 in 2D problem
-      part_temp%r(3) = part%r(3) + alpha
       IF (r > 0.D0) THEN
         sin_alpha = y_new/r
         cos_alpha = x_new/r
@@ -57,7 +95,47 @@ MODULE moduleSolver
       !Copy temporal particle to particle
       part=part_temp
 
-    END SUBROUTINE push
+    END SUBROUTINE pushCylNeutral
+
+    !Push one particle. Boris pusher for 2D Cyl Charged particle
+    PURE SUBROUTINE pushCylCharged(part)
+      USE moduleSpecies
+      USE moduleEM
+      IMPLICIT NONE
+
+      TYPE(particle), INTENT(inout):: part
+      REAL(8):: v_p_oh_star(2:3)
+      TYPE(particle):: part_temp
+      REAL(8):: x_new, y_new, r, sin_alpha, cos_alpha
+      REAL(8):: qmEFt(1:3)!charge*tau*EF/mass
+
+      part_temp = part
+      !Get electric field at particle position
+      qmEFt = part_temp%qm*gatherElecField(part_temp)*tau
+      !z
+      part_temp%v(1) = part%v(1) + qmEFt(1)
+      part_temp%r(1) = part%r(1) + part_temp%v(1)*tau
+      !r,theta
+      v_p_oh_star(2) = part%v(2) + qmEFt(2)
+      x_new = part%r(2) + v_p_oh_star(2)*tau
+      v_p_oh_star(3) = part%v(3) + qmEFt(3)
+      y_new =             v_p_oh_star(3)*tau
+      r = DSQRT(x_new**2+y_new**2)
+      part_temp%r(2) = r
+      IF (r > 0.D0) THEN
+        sin_alpha = y_new/r
+        cos_alpha = x_new/r
+      ELSE
+        sin_alpha = 0.D0
+        cos_alpha = 1.D0
+      END IF
+      part_temp%v(2) =  cos_alpha*v_p_oh_star(2)+sin_alpha*v_p_oh_star(3)
+      part_temp%v(3) = -sin_alpha*v_p_oh_star(2)+cos_alpha*v_p_oh_star(3)
+      part_temp%n_in = .FALSE. !Assume particle is outside until cell is found
+      !Copy temporal particle to particle
+      part=part_temp
+
+    END SUBROUTINE pushCylCharged
 
     !Do the collisions in all the cells
     SUBROUTINE doCollisions()
@@ -153,6 +231,64 @@ MODULE moduleSolver
 
     END SUBROUTINE doScatter
 
+    SUBROUTINE doEMField()
+      IMPLICIT NONE
+
+      CALL solver%solveEMField()
+
+    END SUBROUTINE doEMField
+
+    !Solving the Poisson equation for electrostatic potential
+    SUBROUTINE solveElecField()
+      USE moduleEM
+      USE moduleMesh
+      USE moduleErrors
+      IMPLICIT NONE
+
+      INTEGER, SAVE:: INFO
+      INTEGER:: n
+      REAL(8), ALLOCATABLE, SAVE:: tempF(:)
+      EXTERNAL:: dgetrs
+
+      !$OMP SINGLE
+      ALLOCATE(tempF(1:mesh%numNodes))
+      !$OMP END SINGLE
+
+      CALL assembleSourceVector(tempF)
+
+      !$OMP SINGLE
+      CALL dgetrs('N', mesh%numNodes, 1, mesh%K, mesh%numNodes, &
+                              mesh%IPIV, tempF,  mesh%numNodes, info)
+      !$OMP END SINGLE
+
+      IF (info == 0) THEN
+        !Suscessful resolution of Poission equation
+        !$OMP DO
+        DO n = 1, mesh%numNodes
+          mesh%nodes(n)%obj%emData%phi = tempF(n)
+
+        END DO
+        !$OMP END DO
+
+      ELSE
+        !$OMP SINGLE
+        CALL criticalError('Poisson equation failed', 'solveElecField')
+        !$OMP END SINGLE
+
+      END IF
+
+      !$OMP SINGLE
+      DEALLOCATE(tempF)
+      !$OMP END SINGLE
+
+    END SUBROUTINE solveElecField
+
+    !Empty module that does no computation of EM field for neutral cases
+    SUBROUTINE noEMField()
+      IMPLICIT NONE
+
+    END SUBROUTINE noEMField
+
     SUBROUTINE doOutput(t)
       USE moduleMesh
       USE moduleOutput
@@ -172,8 +308,9 @@ MODULE moduleSolver
         CALL mesh%printOutput(t)
         CALL printTime(t, t == 0)
         CALL mesh%printColl(t)
-        WRITE(*, "(5X,A21,I10,A1,I8)") "t/tmax: ",    t,       "/", tmax
-        WRITE(*, "(5X,A21,I10)")       "Particles: ", nPartOld
+        CALL mesh%printEM(t)
+        WRITE(*, "(5X,A21,I10,A1,I10)") "t/tmax: ",    t,       "/", tmax
+        WRITE(*, "(5X,A21,I10)")        "Particles: ", nPartOld
         IF (t == 0) THEN
           WRITE(*, "(5X,A21,F8.1,A2)")   "     init time: ",    1.D3*tStep, "ms"
 
@@ -183,7 +320,7 @@ MODULE moduleSolver
         END IF
 
         IF (nPartOld > 0) THEN
-          WRITE(*, "(5X,A21,F8.1,A2)") "avg t/particle: ", 1.D9*(tPush+tReset+tColl+tWeight)/DBLE(nPartOld), "ns"
+          WRITE(*, "(5X,A21,F8.1,A2)") "avg t/particle: ", 1.D9*tStep/DBLE(nPartOld), "ns"
 
         END IF
         WRITE(*,*)
diff --git a/src/modules/moduleSpecies.f95 b/src/modules/moduleSpecies.f95
index 9cf9875..dce36a4 100644
--- a/src/modules/moduleSpecies.f95
+++ b/src/modules/moduleSpecies.f95
@@ -31,9 +31,10 @@ MODULE moduleSpecies
     REAL(8):: v(1:3) !Velocity
     INTEGER:: pt !Particle species id
     INTEGER:: e_p !Index of element in which the particle is located
-    REAL(8):: xLog(1:2) !Logical coordinates of particle in element e_p.
+    REAL(8):: xLog(1:3) !Logical coordinates of particle in element e_p.
     LOGICAL:: n_in !Flag that indicates if a particle is in the domain
-    REAL(8):: weight=0.D0
+    REAL(8):: weight=0.D0 !weight of particle
+    REAL(8):: qm = 0.D0 !charge over mass
 
   END TYPE particle