diff --git a/src/main.cpp b/src/main.cpp
index b6d83c7..2930784 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -12,143 +12,84 @@
 
 #include <fstream>
 #include <omp.h>
+#include <string>
 #include <sys/stat.h>
+#include <vector>
 
 #include "PenningTrap.hpp"
+#include "utils.hpp"
 
 #define PARTICLES 100
 #define N 10000
 #define CHARGE 1.
 #define MASS 40. // unit: amu
 
-Particle p1(1., 40., arma::vec{20., 0., 20.}, arma::vec{0., 25., 0.});
-Particle p2(1., 40., arma::vec{25., 25., 0.}, arma::vec{0., 40., 5.});
+Particle p1(CHARGE, MASS, arma::vec{20., 0., 20.}, arma::vec{0., 25., 0.});
+Particle p2(CHARGE, MASS, arma::vec{25., 25., 0.}, arma::vec{0., 40., 5.});
 
 void simulate_single_particle()
 {
-    PenningTrap trap;
-
-    trap.add_particle(p1);
+    PenningTrap trap(std::vector<Particle>{p1});
 
     double time = 50.; // microseconds
-    double dt = time / (double)N;
-    
-    auto res = new arma::vec::fixed<3>[N];
 
-    for (int i=0; i<N;i++) {
-        res[i] = trap.get_particle(0);
-        trap.evolve_RK4(dt);
-    }
-
-    mkdir("output", 0777);
-    mkdir("output/simulate_single_particle", 0777);
-
-    std::ofstream ofile;
-
-    ofile.open("output/simulate_single_particle/out.txt");
-    for (int i=0; i<N;i++) {
-        arma::vec p = res[i];
-        ofile << p(0) << "," << p(1) << "," << p(3) << "\n";
-    }
-    ofile.close();
-    delete [] res;
+    trap.write_simulation_to_dir("output/simulate_single_particle", time, N);
 }
 
 void simulate_two_particles()
 {
-    PenningTrap trap;
-
-    trap.add_particle(p1);
-    trap.add_particle(p2);
+    PenningTrap trap_no_interaction(std::vector<Particle>{p1, p2});
+    PenningTrap trap_with_interaction(std::vector<Particle>{p1, p2});
 
     double time = 50.; // microseconds
-    double dt = time / (double)N;
 
-    auto res = new arma::vec::fixed<3>[2][N];
+    trap_no_interaction.write_simulation_to_dir(
+        "output/simulate_2_particles/no_interaction", time, N, "rk4", false);
+    trap_with_interaction.write_simulation_to_dir(
+        "output/simulate_2_particles/with_interaction", time, N);
+}
 
+void simulate_single_particle_with_different_steps()
+{
 
-    int counter = 0;
+    double time = 50; // microseconds
 
-    // Get the path of all particles
-    for (int j = 0; j < N; j++) {
-#pragma omp parallel for
-        for (int i = 0; i < 2; i++) {
-            res[i][j] = trap.get_particle(i);
-        }
-        trap.evolve_RK4(dt);
+    for (int i = 0; i < 4; i++) {
+        int steps = 4000 * (i + 1);
+        PenningTrap trap(std::vector<Particle>{p1});
+        trap.write_simulation_to_dir("output/N_steps/RK4/" + std::to_string(steps) +
+                                         "_steps",
+                                     time, steps, "rk4", false);
     }
 
-    std::cout << counter << std::endl;
-
-    arma::vec::fixed<3> *cur_row;
-    arma::vec::fixed<3> cur_elem;
-
-    mkdir("output", 0777);
-    mkdir("output/simulate_2_particles", 0777);
-
-    std::ofstream ofile;
-
-// Write particle paths to file
-#pragma omp parallel for private(cur_row, cur_elem, ofile)
-    for (int i = 0; i < 2; i++) {
-        cur_row = res[i];
-        ofile.open("output/simulate_100_particles/p" + std::to_string(i) +
-                   ".txt");
-        for (int j = 0; j < N; j++) {
-            cur_elem = cur_row[j];
-            ofile << cur_elem(0) << "," << cur_elem(1) << "," << cur_elem(2)
-                  << "\n";
-        }
-        ofile.close();
+    for (int i = 0; i < 4; i++) {
+        int steps = 4000 * (i + 1);
+        PenningTrap trap(std::vector<Particle>{p1});
+        trap.write_simulation_to_dir("output/N_steps/euler/" + std::to_string(steps) +
+                                         "_steps",
+                                     time, steps, "euler", false);
     }
-
 }
 
 void simulate_100_particles()
 {
-    PenningTrap trap;
-
-    // Add particles inside trap
-    for (int i = 0; i < PARTICLES; i++) {
-        arma::vec r = arma::vec(3).randn() * 0.1 *
-                      trap.get_d(); // random initial position
-        arma::vec v = arma::vec(3).randn() * 0.1 *
-                      trap.get_d(); // random initial velocity
-        trap.add_particle(Particle(CHARGE, MASS, r, v));
-    }
+    PenningTrap trap(100);
 
     double time = 50.; // microseconds
-    double dt = time / (double)N;
 
-    sim_arr res = trap.simulate(time, N);
-
-    sim_rows cur_row;
-    arma::vec::fixed<3> cur_elem;
-
-    mkdir("output", 0777);
-    mkdir("output/simulate_100_particles", 0777);
-
-    std::ofstream ofile;
-
-// Write particle paths to file
-#pragma omp parallel for private(cur_row, cur_elem, ofile)
-    for (int i = 0; i < PARTICLES; i++) {
-        cur_row = res[i];
-        ofile.open("output/simulate_100_particles/p" + std::to_string(i) +
-                   ".txt");
-        for (int j = 0; j < N; j++) {
-            cur_elem = cur_row[j];
-            ofile << cur_elem(0) << "," << cur_elem(1) << "," << cur_elem(2)
-                  << "\n";
-        }
-        ofile.close();
-    }
+    trap.write_simulation_to_dir("output/simulate_100_particles", time, N);
 }
 
 int main()
 {
     double start = omp_get_wtime();
 
+    simulate_single_particle();
+
+    simulate_two_particles();
+
+    simulate_single_particle_with_different_steps();
+
     simulate_100_particles();
 
     double end = omp_get_wtime();