diff --git a/src/makefile b/src/makefile
index 24ad089..2a421f8 100644
--- a/src/makefile
+++ b/src/makefile
@@ -1,4 +1,4 @@
-OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleMeshBoundary.o $(OBJDIR)/moduleCompTime.o $(OBJDIR)/moduleSolver.o \
+OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleMeshBoundary.o $(OBJDIR)/moduleCompTime.o \
 			  	$(OBJDIR)/moduleSpecies.o $(OBJDIR)/moduleInject.o $(OBJDIR)/moduleInput.o \
 			  	$(OBJDIR)/moduleErrors.o $(OBJDIR)/moduleList.o $(OBJDIR)/moduleOutput.o \
 			  	$(OBJDIR)/moduleBoundary.o $(OBJDIR)/moduleCaseParam.o $(OBJDIR)/moduleRefParam.o \
@@ -12,7 +12,9 @@ OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleMeshBoundary.o $(OBJDIR)/module
 			  	$(OBJDIR)/moduleMesh2DCart.o \
           $(OBJDIR)/moduleMesh1DRad.o  \
           $(OBJDIR)/moduleMesh1DCart.o \
-          $(OBJDIR)/moduleMesh0D.o
+          $(OBJDIR)/moduleMesh0D.o     \
+          $(OBJDIR)/moduleSolver.o     \
+					$(OBJDIR)/modulePusher.o
 
 
 all: $(OUTPUT)
diff --git a/src/modules/common/makefile b/src/modules/common/makefile
new file mode 100644
index 0000000..a57871b
--- /dev/null
+++ b/src/modules/common/makefile
@@ -0,0 +1,11 @@
+OBJS = moduleCompTime.o moduleCaseParam.o moduleConstParam.o \
+			 moduleErrors.o moduleMath.o moduleParallel.o \
+			 moduleRandom.o moduleRefParam.o moduleTable.o
+
+all: $(OBJS)
+
+moduleTable.o: moduleErrors.o moduleTable.f90
+	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
+
+%.o: %.f90
+	$(FC) $(FCFLAGS) -c $< -o $(OBJDIR)/$@
diff --git a/src/modules/moduleCaseParam.f90 b/src/modules/common/moduleCaseParam.f90
similarity index 100%
rename from src/modules/moduleCaseParam.f90
rename to src/modules/common/moduleCaseParam.f90
diff --git a/src/modules/moduleCompTime.f90 b/src/modules/common/moduleCompTime.f90
similarity index 91%
rename from src/modules/moduleCompTime.f90
rename to src/modules/common/moduleCompTime.f90
index d0a6755..46310f0 100644
--- a/src/modules/moduleCompTime.f90
+++ b/src/modules/common/moduleCompTime.f90
@@ -1,5 +1,9 @@
 !Information to calculate computation time
 MODULE moduleCompTime
+  IMPLICIT NONE
+
+  PUBLIC
+
   REAL(8):: tStep    = 0.D0
   REAL(8):: tPush    = 0.D0
   REAL(8):: tReset   = 0.D0
diff --git a/src/modules/moduleConstParam.f90 b/src/modules/common/moduleConstParam.f90
similarity index 100%
rename from src/modules/moduleConstParam.f90
rename to src/modules/common/moduleConstParam.f90
diff --git a/src/modules/moduleErrors.f90 b/src/modules/common/moduleErrors.f90
similarity index 100%
rename from src/modules/moduleErrors.f90
rename to src/modules/common/moduleErrors.f90
diff --git a/src/modules/moduleMath.f90 b/src/modules/common/moduleMath.f90
similarity index 100%
rename from src/modules/moduleMath.f90
rename to src/modules/common/moduleMath.f90
diff --git a/src/modules/moduleParallel.f90 b/src/modules/common/moduleParallel.f90
similarity index 100%
rename from src/modules/moduleParallel.f90
rename to src/modules/common/moduleParallel.f90
diff --git a/src/modules/moduleRandom.f90 b/src/modules/common/moduleRandom.f90
similarity index 100%
rename from src/modules/moduleRandom.f90
rename to src/modules/common/moduleRandom.f90
diff --git a/src/modules/moduleRefParam.f90 b/src/modules/common/moduleRefParam.f90
similarity index 100%
rename from src/modules/moduleRefParam.f90
rename to src/modules/common/moduleRefParam.f90
diff --git a/src/modules/moduleTable.f90 b/src/modules/common/moduleTable.f90
similarity index 100%
rename from src/modules/moduleTable.f90
rename to src/modules/common/moduleTable.f90
diff --git a/src/modules/init/makefile b/src/modules/init/makefile
new file mode 100644
index 0000000..c035f7e
--- /dev/null
+++ b/src/modules/init/makefile
@@ -0,0 +1,5 @@
+all: moduleInput.o
+
+%.o: %.f90
+	$(FC) $(FCFLAGS) -c $< -o $(OBJDIR)/$@
+
diff --git a/src/modules/moduleInput.f90 b/src/modules/init/moduleInput.f90
similarity index 99%
rename from src/modules/moduleInput.f90
rename to src/modules/init/moduleInput.f90
index 348edf3..2523b3b 100644
--- a/src/modules/moduleInput.f90
+++ b/src/modules/init/moduleInput.f90
@@ -9,7 +9,6 @@ MODULE moduleInput
       USE json_module
       USE moduleErrors
       USE moduleBoundary
-      USE moduleInject
       USE moduleOutput
       USE moduleMesh
       IMPLICIT NONE
diff --git a/src/modules/makefile b/src/modules/makefile
index f1545dc..254d015 100644
--- a/src/modules/makefile
+++ b/src/modules/makefile
@@ -1,38 +1,40 @@
-
-OBJS = moduleCaseParam.o moduleCompTime.o moduleList.o \
-		   moduleOutput.o moduleInput.o moduleSolver.o \
-			 moduleCollisions.o moduleTable.o moduleParallel.o \
-			 moduleEM.o moduleRandom.o moduleMath.o \
-			 moduleProbe.o
+OBJS = common.o output.o mesh.o solver.o init.o \
+			 moduleBoundary.o moduleCollisions.o moduleInject.o \
+			 moduleList.o moduleProbe.o \
+			 moduleSpecies.o
 
 all: $(OBJS)
 
-mesh.o: moduleCollisions.o moduleBoundary.o moduleAverage.o
+common.o:
+	$(MAKE) -C common all
+
+output.o: moduleSpecies.o common.o
+	$(MAKE) -C output all
+
+mesh.o: moduleCollisions.o moduleBoundary.o output.o
 	$(MAKE) -C mesh all
 
-moduleCollisions.o: moduleList.o moduleMath.o moduleRandom.o moduleTable.o moduleSpecies.o moduleRefParam.o moduleConstParam.o moduleCollisions.f90
+solver.o: moduleSpecies.o moduleProbe.o common.o output.o mesh.o
+	$(MAKE) -C solver all
+
+init.o: common.o solver.o moduleInject.o
+	$(MAKE) -C init all
+
+moduleBoundary.o: common.o moduleBoundary.f90
 	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
 
-moduleInput.o: moduleParallel.o moduleSolver.o moduleInject.o moduleInput.f90
+moduleCollisions.o: moduleList.o moduleSpecies.o common.o moduleCollisions.f90
 	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
 
-moduleInject.o: moduleInject.f90
-	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
-
-moduleList.o: moduleConstParam.o moduleErrors.o moduleCaseParam.o moduleSpecies.o moduleList.f90
-	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
-
-moduleOutput.o: moduleMath.o moduleRefParam.o moduleOutput.f90
-	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
-
-moduleSolver.o: moduleProbe.o moduleEM.o moduleAverage.o moduleSolver.f90
+moduleList.o: common.o moduleSpecies.o moduleList.f90
 	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
 
 moduleProbe.o: mesh.o moduleProbe.f90
 	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
 
-moduleBoundary.o: moduleTable.o moduleBoundary.f90
+moduleSpecies.o: common.o moduleSpecies.f90
 	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
 
 %.o: %.f90
 	$(FC) $(FCFLAGS) -c $< -o $(OBJDIR)/$@
+
diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90
index 660a313..b8edfdd 100644
--- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90
+++ b/src/modules/mesh/1DRad/moduleMesh1DRad.f90
@@ -265,11 +265,11 @@ MODULE moduleMesh1DRad
       l = 2.D0*detJ
       self%volume  = r*l
       !Computes volume per node
-      xi = (/-5.D-1,  0.D0, 0.D0/)
-      r = DOT_PRODUCT(self%fPsi(xi),self%r)
+      Xii = (/-5.D-1,  0.D0, 0.D0/)
+      r = DOT_PRODUCT(self%fPsi(Xii),self%r)
       self%arNodes(1) = fPsi(1)*r*l
-      xi = (/ 5.D-1,  0.D0, 0.D0/)
-      r = DOT_PRODUCT(self%fPsi(xi),self%r)
+      Xii = (/ 5.D-1,  0.D0, 0.D0/)
+      r = DOT_PRODUCT(self%fPsi(Xii),self%r)
       self%arNodes(2) = fPsi(2)*r*l
 
     END SUBROUTINE areaRad
diff --git a/src/modules/mesh/makefile b/src/modules/mesh/makefile
index a013878..7c66c3b 100644
--- a/src/modules/mesh/makefile
+++ b/src/modules/mesh/makefile
@@ -26,3 +26,6 @@ moduleMeshBoundary.o: moduleMesh.o moduleMeshBoundary.f90
 
 inout.o: 3DCart.o 2DCyl.o 2DCart.o 1DRad.o 1DCart.o 0D.o
 	$(MAKE) -C inout all
+
+%.o: %.f90
+	$(FC) $(FCFLAGS) -c $< -o $(OBJDIR)/$@
diff --git a/src/modules/output/makefile b/src/modules/output/makefile
new file mode 100644
index 0000000..df42742
--- /dev/null
+++ b/src/modules/output/makefile
@@ -0,0 +1,7 @@
+all: moduleAverage.o moduleOutput.o
+
+moduleAverage.o: moduleOutput.o moduleAverage.f90
+	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
+
+%.o: %.f90
+	$(FC) $(FCFLAGS) -c $< -o $(OBJDIR)/$@
diff --git a/src/modules/moduleAverage.f90 b/src/modules/output/moduleAverage.f90
similarity index 100%
rename from src/modules/moduleAverage.f90
rename to src/modules/output/moduleAverage.f90
diff --git a/src/modules/moduleOutput.f90 b/src/modules/output/moduleOutput.f90
similarity index 100%
rename from src/modules/moduleOutput.f90
rename to src/modules/output/moduleOutput.f90
diff --git a/src/modules/solver/electromagnetic/makefile b/src/modules/solver/electromagnetic/makefile
new file mode 100644
index 0000000..e4dd2fa
--- /dev/null
+++ b/src/modules/solver/electromagnetic/makefile
@@ -0,0 +1,8 @@
+all: moduleEM.o
+
+moduleEM.o: moduleEM.f90
+	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
+
+%.o: %.f90
+	$(FC) $(FCFLAGS) -c $< -o $(OBJDIR)/$@
+
diff --git a/src/modules/moduleEM.f90 b/src/modules/solver/electromagnetic/moduleEM.f90
similarity index 73%
rename from src/modules/moduleEM.f90
rename to src/modules/solver/electromagnetic/moduleEM.f90
index ce91858..ae73ebc 100644
--- a/src/modules/moduleEM.f90
+++ b/src/modules/solver/electromagnetic/moduleEM.f90
@@ -19,6 +19,7 @@ MODULE moduleEM
   REAL(8), ALLOCATABLE:: qSpecies(:)
 
   CONTAINS
+    !Apply boundary conditions to the K matrix for Poisson's equation
     SUBROUTINE apply(self, edge)
       USE moduleMesh
       IMPLICIT NONE
@@ -79,6 +80,7 @@ MODULE moduleEM
 
     END FUNCTION gatherMagnField
 
+    !Assemble the source vector based on the charge density to solve Poisson's equation
     SUBROUTINE assembleSourceVector(vectorF)
       USE moduleMesh
       USE moduleRefParam
@@ -142,4 +144,48 @@ MODULE moduleEM
 
     END SUBROUTINE assembleSourceVector
 
+    !Solving the Poisson equation for electrostatic potential
+    SUBROUTINE solveElecField()
+      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
+
 END MODULE moduleEM
diff --git a/src/modules/solver/makefile b/src/modules/solver/makefile
new file mode 100644
index 0000000..d50bdc8
--- /dev/null
+++ b/src/modules/solver/makefile
@@ -0,0 +1,13 @@
+all: moduleSolver.o electromagnetic.o pusher.o
+
+electromagnetic.o:
+	$(MAKE) -C electromagnetic all
+
+pusher.o:
+	$(MAKE) -C pusher all
+
+moduleSolver.o: electromagnetic.o pusher.o moduleSolver.f90
+	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
+
+%.o: %.f90
+	$(FC) $(FCFLAGS) -c $< -o $(OBJDIR)/$@
diff --git a/src/modules/moduleSolver.f90 b/src/modules/solver/moduleSolver.f90
similarity index 64%
rename from src/modules/moduleSolver.f90
rename to src/modules/solver/moduleSolver.f90
index d0b2d01..58a7a71 100644
--- a/src/modules/moduleSolver.f90
+++ b/src/modules/solver/moduleSolver.f90
@@ -62,6 +62,7 @@ MODULE moduleSolver
     !Init Pusher
     SUBROUTINE initPusher(self, pusherType, tau, tauSp)
       USE moduleErrors
+      USE modulePusher
       USE moduleMesh, ONLY: mesh
       IMPLICIT NONE
 
@@ -128,6 +129,7 @@ MODULE moduleSolver
     END SUBROUTINE initPusher
 
     SUBROUTINE initEM(self, EMType)
+      USE moduleEM
       IMPLICIT NONE
 
       CLASS(solverGeneric), INTENT(inout):: self
@@ -184,239 +186,6 @@ MODULE moduleSolver
 
     END SUBROUTINE doPushes
 
-    !Push neutral particles in cartesian coordinates
-    PURE SUBROUTINE pushCartNeutral(part, tauIn)
-      USE moduleSPecies
-      IMPLICIT NONE
-
-      TYPE(particle), INTENT(inout):: part
-      REAL(8), INTENT(in):: tauIn
-
-      part%r = part%r + part%v*tauIn
-
-    END SUBROUTINE pushCartNeutral
-
-    PURE SUBROUTINE pushCartElectrostatic(part, tauIn)
-      USE moduleSPecies
-      USE moduleEM
-      IMPLICIT NONE
-
-      TYPE(particle), INTENT(inout):: part
-      REAL(8), INTENT(in):: tauIn
-      REAL(8):: qmEFt(1:3)
-
-      !Get the electric field at particle position
-      qmEFt = part%species%qm*gatherElecField(part)*tauIn
-
-      !Update velocity
-      part%v = part%v + qmEFt
-
-      !Update position
-      part%r = part%r + part%v*tauIn
-
-    END SUBROUTINE pushCartElectrostatic
-
-    PURE SUBROUTINE pushCartElectromagnetic(part, tauIn)
-      USE moduleSPecies
-      USE moduleEM
-      USE moduleMath
-      IMPLICIT NONE
-
-      TYPE(particle), INTENT(inout):: part
-      REAL(8), INTENT(in):: tauIn
-      REAL(8):: tauInHalf
-      REAL(8):: qmEFt(1:3)
-      REAL(8):: B(1:3), BNorm
-      REAL(8):: fn
-      REAL(8):: v_minus(1:3), v_prime(1:3), v_plus(1:3)
-
-      tauInHalf = tauIn *0.5D0
-      !Half of the force o f the electric field
-      qmEFt = part%species%qm*gatherElecField(part)*tauInHalf
-
-      !Half step for electrostatic
-      v_minus = part%v + qmEFt
-
-      !Full step rotation
-      B = gatherMagnField(part)
-      BNorm = NORM2(B)
-      IF (BNorm > 0.D0) THEN
-        fn = DTAN(part%species%qm * tauInHalf*BNorm) / BNorm
-        v_prime = v_minus + fn * crossProduct(v_minus, B)
-        v_plus  = v_minus + 2.D0 * fn / (1.D0 + fn**2 * B**2)*crossProduct(v_prime, B)
-
-      END IF
-
-      !Half step for electrostatic
-      part%v = v_plus + qmEFt
-
-      !Update position
-      part%r = part%r + part%v*tauIn
-
-    END SUBROUTINE pushCartElectromagnetic
-
-    !Push one particle. Boris pusher for 2D Cyl Neutral particle
-    PURE SUBROUTINE push2DCylNeutral(part, tauIn)
-      USE moduleSpecies
-      IMPLICIT NONE
-
-      TYPE(particle), INTENT(inout):: part
-      REAL(8), INTENT(in):: tauIn
-      TYPE(particle):: part_temp
-      REAL(8):: x_new, y_new, r, sin_alpha, cos_alpha
-      REAL(8):: v_p_oh_star(2:3)
-
-      part_temp = part
-      !z
-      part_temp%v(1) = part%v(1)
-      part_temp%r(1) = part%r(1) + part_temp%v(1)*tauIn
-      !r,theta
-      v_p_oh_star(2) = part%v(2)
-      x_new = part%r(2) + v_p_oh_star(2)*tauIn
-      v_p_oh_star(3) = part%v(3)
-      y_new =             v_p_oh_star(3)*tauIn
-      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)
-
-      !Copy temporal particle to particle
-      part=part_temp
-
-    END SUBROUTINE push2DCylNeutral
-
-    !Push one particle. Boris pusher for 2D Cyl Electrostatic particle
-    PURE SUBROUTINE push2DCylElectrostatic(part, tauIn)
-      USE moduleSpecies
-      USE moduleEM
-      IMPLICIT NONE
-
-      TYPE(particle), INTENT(inout):: part
-      REAL(8), INTENT(in):: tauIn
-      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*tauIn*EF/mass
-
-      part_temp = part
-      !Get electric field at particle position
-      qmEFt = part_temp%species%qm*gatherElecField(part_temp)*tauIn
-      !z
-      part_temp%v(1) = part%v(1) + qmEFt(1)
-      part_temp%r(1) = part%r(1) + part_temp%v(1)*tauIn
-      !r,theta
-      v_p_oh_star(2) = part%v(2) + qmEFt(2)
-      x_new = part%r(2) + v_p_oh_star(2)*tauIn
-      v_p_oh_star(3) = part%v(3) + qmEFt(3)
-      y_new =             v_p_oh_star(3)*tauIn
-      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)
-
-      !Copy temporal particle to particle
-      part=part_temp
-
-    END SUBROUTINE push2DCylElectrostatic
-
-    !Push one particle. Boris pusher for 1D Radial Neutral particle
-    PURE SUBROUTINE push1DRadNeutral(part, tauIn)
-      USE moduleSpecies
-      USE moduleEM
-      IMPLICIT NONE
-
-      TYPE(particle), INTENT(inout):: part
-      REAL(8), INTENT(in):: tauIn
-      REAL(8):: v_p_oh_star(1:2)
-      TYPE(particle):: part_temp
-      REAL(8):: x_new, y_new, r, sin_alpha, cos_alpha
-
-      part_temp = part
-      !r,theta
-      v_p_oh_star(1) = part%v(1)
-      x_new = part%r(1) + v_p_oh_star(1)*tauIn
-      v_p_oh_star(2) = part%v(2)
-      y_new =             v_p_oh_star(2)*tauIn
-      r = DSQRT(x_new**2+y_new**2)
-      part_temp%r(1) = 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(1) =  cos_alpha*v_p_oh_star(1)+sin_alpha*v_p_oh_star(2)
-      part_temp%v(2) = -sin_alpha*v_p_oh_star(1)+cos_alpha*v_p_oh_star(2)
-
-      !Copy temporal particle to particle
-      part=part_temp
-
-    END SUBROUTINE push1DRadNeutral
-
-    !Push one particle. Boris pusher for 1D Radial Electrostatic particle
-    PURE SUBROUTINE push1DRadElectrostatic(part, tauIn)
-      USE moduleSpecies
-      USE moduleEM
-      IMPLICIT NONE
-
-      TYPE(particle), INTENT(inout):: part
-      REAL(8), INTENT(in):: tauIn
-      REAL(8):: v_p_oh_star(1:2)
-      TYPE(particle):: part_temp
-      REAL(8):: x_new, y_new, r, sin_alpha, cos_alpha
-      REAL(8):: qmEFt(1:3)!charge*tauIn*EF/mass
-
-      part_temp = part
-      !Get electric field at particle position
-      qmEFt = part_temp%species%qm*gatherElecField(part_temp)*tauMin
-      !r,theta
-      v_p_oh_star(1) = part%v(1) + qmEFt(1)
-      x_new = part%r(1) + v_p_oh_star(1)*tauIn
-      v_p_oh_star(2) = part%v(2) + qmEFt(2)
-      y_new =             v_p_oh_star(2)*tauIn
-      r = DSQRT(x_new**2+y_new**2)
-      part_temp%r(1) = 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(1) =  cos_alpha*v_p_oh_star(1)+sin_alpha*v_p_oh_star(2)
-      part_temp%v(2) = -sin_alpha*v_p_oh_star(1)+cos_alpha*v_p_oh_star(2)
-
-      !Copy temporal particle to particle
-      part=part_temp
-
-    END SUBROUTINE push1DRadElectrostatic
-
-    !Dummy pusher for 0D geometry
-    PURE SUBROUTINE push0D(part, tauIn)
-      USE moduleSpecies
-      USE moduleEM
-      IMPLICIT NONE
-
-      TYPE(particle), INTENT(inout):: part
-      REAL(8), INTENT(in):: tauIn
-
-    END SUBROUTINE push0D
-
     !Takes the particles from a list and put them into an array
     !nStart indicates the last fill index in the array
     SUBROUTINE resetList(partList, partArray, nStart)
@@ -606,55 +375,11 @@ MODULE moduleSolver
 
     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
-
-    !Modify particle weight as a function of cell volume and splits particle
+    !Split particles as a function of cell volume and splits particle
     SUBROUTINE volumeWScheme(part, volOld, volNew)
       USE moduleSpecies
       USE moduleMesh
+      USE moduleRandom
       IMPLICIT NONE
 
       TYPE(particle), INTENT(inout):: part
@@ -669,7 +394,7 @@ MODULE moduleSolver
         !Calculate probability of splitting particle
         pSplit = 1.D0 - DEXP(-fractionVolume)
 
-        IF (random() < pSplit THEN
+        IF (random() < pSplit) THEN
           !Split particle in two
           CALL splitParticle(part, 2, volNew)
 
diff --git a/src/modules/solver/pusher/makefile b/src/modules/solver/pusher/makefile
new file mode 100644
index 0000000..bf574d5
--- /dev/null
+++ b/src/modules/solver/pusher/makefile
@@ -0,0 +1,4 @@
+all: modulePusher.o
+
+%.o: %.f90
+	$(FC) $(FCFLAGS) -c $< -o $(OBJDIR)/$@
diff --git a/src/modules/solver/pusher/modulePusher.f90 b/src/modules/solver/pusher/modulePusher.f90
new file mode 100644
index 0000000..69fcaab
--- /dev/null
+++ b/src/modules/solver/pusher/modulePusher.f90
@@ -0,0 +1,237 @@
+MODULE modulePusher
+
+  CONTAINS
+    !Push neutral particles in cartesian coordinates
+    PURE SUBROUTINE pushCartNeutral(part, tauIn)
+      USE moduleSpecies
+      IMPLICIT NONE
+
+      TYPE(particle), INTENT(inout):: part
+      REAL(8), INTENT(in):: tauIn
+
+      part%r = part%r + part%v*tauIn
+
+    END SUBROUTINE pushCartNeutral
+
+    PURE SUBROUTINE pushCartElectrostatic(part, tauIn)
+      USE moduleSPecies
+      USE moduleEM
+      IMPLICIT NONE
+
+      TYPE(particle), INTENT(inout):: part
+      REAL(8), INTENT(in):: tauIn
+      REAL(8):: qmEFt(1:3)
+
+      !Get the electric field at particle position
+      qmEFt = part%species%qm*gatherElecField(part)*tauIn
+
+      !Update velocity
+      part%v = part%v + qmEFt
+
+      !Update position
+      part%r = part%r + part%v*tauIn
+
+    END SUBROUTINE pushCartElectrostatic
+
+    PURE SUBROUTINE pushCartElectromagnetic(part, tauIn)
+      USE moduleSPecies
+      USE moduleEM
+      USE moduleMath
+      IMPLICIT NONE
+
+      TYPE(particle), INTENT(inout):: part
+      REAL(8), INTENT(in):: tauIn
+      REAL(8):: tauInHalf
+      REAL(8):: qmEFt(1:3)
+      REAL(8):: B(1:3), BNorm
+      REAL(8):: fn
+      REAL(8):: v_minus(1:3), v_prime(1:3), v_plus(1:3)
+
+      tauInHalf = tauIn *0.5D0
+      !Half of the force o f the electric field
+      qmEFt = part%species%qm*gatherElecField(part)*tauInHalf
+
+      !Half step for electrostatic
+      v_minus = part%v + qmEFt
+
+      !Full step rotation
+      B = gatherMagnField(part)
+      BNorm = NORM2(B)
+      IF (BNorm > 0.D0) THEN
+        fn = DTAN(part%species%qm * tauInHalf*BNorm) / BNorm
+        v_prime = v_minus + fn * crossProduct(v_minus, B)
+        v_plus  = v_minus + 2.D0 * fn / (1.D0 + fn**2 * B**2)*crossProduct(v_prime, B)
+
+      END IF
+
+      !Half step for electrostatic
+      part%v = v_plus + qmEFt
+
+      !Update position
+      part%r = part%r + part%v*tauIn
+
+    END SUBROUTINE pushCartElectromagnetic
+
+    !Push one particle. Boris pusher for 2D Cyl Neutral particle
+    PURE SUBROUTINE push2DCylNeutral(part, tauIn)
+      USE moduleSpecies
+      IMPLICIT NONE
+
+      TYPE(particle), INTENT(inout):: part
+      REAL(8), INTENT(in):: tauIn
+      TYPE(particle):: part_temp
+      REAL(8):: x_new, y_new, r, sin_alpha, cos_alpha
+      REAL(8):: v_p_oh_star(2:3)
+
+      part_temp = part
+      !z
+      part_temp%v(1) = part%v(1)
+      part_temp%r(1) = part%r(1) + part_temp%v(1)*tauIn
+      !r,theta
+      v_p_oh_star(2) = part%v(2)
+      x_new = part%r(2) + v_p_oh_star(2)*tauIn
+      v_p_oh_star(3) = part%v(3)
+      y_new =             v_p_oh_star(3)*tauIn
+      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)
+
+      !Copy temporal particle to particle
+      part=part_temp
+
+    END SUBROUTINE push2DCylNeutral
+
+    !Push one particle. Boris pusher for 2D Cyl Electrostatic particle
+    PURE SUBROUTINE push2DCylElectrostatic(part, tauIn)
+      USE moduleSpecies
+      USE moduleEM
+      IMPLICIT NONE
+
+      TYPE(particle), INTENT(inout):: part
+      REAL(8), INTENT(in):: tauIn
+      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*tauIn*EF/mass
+
+      part_temp = part
+      !Get electric field at particle position
+      qmEFt = part_temp%species%qm*gatherElecField(part_temp)*tauIn
+      !z
+      part_temp%v(1) = part%v(1) + qmEFt(1)
+      part_temp%r(1) = part%r(1) + part_temp%v(1)*tauIn
+      !r,theta
+      v_p_oh_star(2) = part%v(2) + qmEFt(2)
+      x_new = part%r(2) + v_p_oh_star(2)*tauIn
+      v_p_oh_star(3) = part%v(3) + qmEFt(3)
+      y_new =             v_p_oh_star(3)*tauIn
+      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)
+
+      !Copy temporal particle to particle
+      part=part_temp
+
+    END SUBROUTINE push2DCylElectrostatic
+
+    !Push one particle. Boris pusher for 1D Radial Neutral particle
+    PURE SUBROUTINE push1DRadNeutral(part, tauIn)
+      USE moduleSpecies
+      USE moduleEM
+      IMPLICIT NONE
+
+      TYPE(particle), INTENT(inout):: part
+      REAL(8), INTENT(in):: tauIn
+      REAL(8):: v_p_oh_star(1:2)
+      TYPE(particle):: part_temp
+      REAL(8):: x_new, y_new, r, sin_alpha, cos_alpha
+
+      part_temp = part
+      !r,theta
+      v_p_oh_star(1) = part%v(1)
+      x_new = part%r(1) + v_p_oh_star(1)*tauIn
+      v_p_oh_star(2) = part%v(2)
+      y_new =             v_p_oh_star(2)*tauIn
+      r = DSQRT(x_new**2+y_new**2)
+      part_temp%r(1) = 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(1) =  cos_alpha*v_p_oh_star(1)+sin_alpha*v_p_oh_star(2)
+      part_temp%v(2) = -sin_alpha*v_p_oh_star(1)+cos_alpha*v_p_oh_star(2)
+
+      !Copy temporal particle to particle
+      part=part_temp
+
+    END SUBROUTINE push1DRadNeutral
+
+    !Push one particle. Boris pusher for 1D Radial Electrostatic particle
+    PURE SUBROUTINE push1DRadElectrostatic(part, tauIn)
+      USE moduleSpecies
+      USE moduleEM
+      IMPLICIT NONE
+
+      TYPE(particle), INTENT(inout):: part
+      REAL(8), INTENT(in):: tauIn
+      REAL(8):: v_p_oh_star(1:2)
+      TYPE(particle):: part_temp
+      REAL(8):: x_new, y_new, r, sin_alpha, cos_alpha
+      REAL(8):: qmEFt(1:3)!charge*tauIn*EF/mass
+
+      part_temp = part
+      !Get electric field at particle position
+      qmEFt = part_temp%species%qm*gatherElecField(part_temp)*tauMin
+      !r,theta
+      v_p_oh_star(1) = part%v(1) + qmEFt(1)
+      x_new = part%r(1) + v_p_oh_star(1)*tauIn
+      v_p_oh_star(2) = part%v(2) + qmEFt(2)
+      y_new =             v_p_oh_star(2)*tauIn
+      r = DSQRT(x_new**2+y_new**2)
+      part_temp%r(1) = 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(1) =  cos_alpha*v_p_oh_star(1)+sin_alpha*v_p_oh_star(2)
+      part_temp%v(2) = -sin_alpha*v_p_oh_star(1)+cos_alpha*v_p_oh_star(2)
+
+      !Copy temporal particle to particle
+      part=part_temp
+
+    END SUBROUTINE push1DRadElectrostatic
+
+    !Dummy pusher for 0D geometry
+    PURE SUBROUTINE push0D(part, tauIn)
+      USE moduleSpecies
+      USE moduleEM
+      IMPLICIT NONE
+
+      TYPE(particle), INTENT(inout):: part
+      REAL(8), INTENT(in):: tauIn
+
+    END SUBROUTINE push0D
+
+END MODULE modulePusher