diff --git a/include/PenningTrap.hpp b/include/PenningTrap.hpp
index eaf89be..9a78784 100644
--- a/include/PenningTrap.hpp
+++ b/include/PenningTrap.hpp
@@ -222,6 +222,9 @@ public:
     double fraction_of_particles_left(double time, unsigned int steps,
                                       std::string method = "rk4",
                                       bool particle_interaction = true);
+
+    vec_3d get_r(int i);
+    double get_t();
 };
 
 #endif
diff --git a/latex/images/3d_plot.pdf b/latex/images/3d_plot.pdf
new file mode 100644
index 0000000..78b115e
Binary files /dev/null and b/latex/images/3d_plot.pdf differ
diff --git a/latex/images/phase_space_2_particles.pdf b/latex/images/phase_space_2_particles.pdf
new file mode 100644
index 0000000..6c006d4
Binary files /dev/null and b/latex/images/phase_space_2_particles.pdf differ
diff --git a/latex/images/phase_space_2_particles_x.pdf b/latex/images/phase_space_2_particles_x.pdf
new file mode 100644
index 0000000..7af3a9c
Binary files /dev/null and b/latex/images/phase_space_2_particles_x.pdf differ
diff --git a/latex/images/phase_space_2_particles_z.pdf b/latex/images/phase_space_2_particles_z.pdf
new file mode 100644
index 0000000..86fe6c3
Binary files /dev/null and b/latex/images/phase_space_2_particles_z.pdf differ
diff --git a/latex/images/plot_2_particles_xy.pdf b/latex/images/plot_2_particles_xy.pdf
new file mode 100644
index 0000000..8e44ccd
Binary files /dev/null and b/latex/images/plot_2_particles_xy.pdf differ
diff --git a/latex/images/single_particle.pdf b/latex/images/single_particle.pdf
index 1180de9..cf0ee5a 100644
Binary files a/latex/images/single_particle.pdf and b/latex/images/single_particle.pdf differ
diff --git a/src/PenningTrap.cpp b/src/PenningTrap.cpp
index ab99599..d5a3b16 100644
--- a/src/PenningTrap.cpp
+++ b/src/PenningTrap.cpp
@@ -53,8 +53,8 @@ vec_3d PenningTrap::v_func(unsigned int i, unsigned int j, double dt)
     case 2:
         return dt * this->k_v[2][j];
     case 3:
-        return (dt / 6.) * (this->k_v[0][j] + 2.*this->k_v[1][j] +
-                            2.*this->k_v[2][j] + this->k_v[3][j]);
+        return (dt / 6.) * (this->k_v[0][j] + 2. * this->k_v[1][j] +
+                            2. * this->k_v[2][j] + this->k_v[3][j]);
     default:
         std::cout << "Not valid!" << std::endl;
         abort();
@@ -71,8 +71,8 @@ vec_3d PenningTrap::r_func(unsigned int i, unsigned int j, double dt)
     case 2:
         return dt * this->k_r[2][j];
     case 3:
-        return (dt / 6.) * (this->k_r[0][j] + 2.*this->k_r[1][j] +
-                            2.*this->k_r[2][j] + this->k_r[3][j]);
+        return (dt / 6.) * (this->k_r[0][j] + 2. * this->k_r[1][j] +
+                            2. * this->k_r[2][j] + this->k_r[3][j]);
     default:
         std::cout << "Not valid!" << std::endl;
         abort();
@@ -106,8 +106,6 @@ vec_3d PenningTrap::force_on_particle(unsigned int i, unsigned int j)
     // Get the distance between the particles
     double norm = arma::norm(res, 2);
 
-    // Multiply res with p_j's charge divided by the norm cubed
-
     return vec_3d(res * p_j.q / (norm * norm * norm));
 }
 
@@ -298,3 +296,13 @@ double PenningTrap::fraction_of_particles_left(double time, unsigned int steps,
 
     return (double)particles_left / (double)this->particles.size();
 }
+
+vec_3d PenningTrap::get_r(int i)
+{
+    return this->particles[i].r_vec;
+}
+
+double PenningTrap::get_t()
+{
+    return this->t;
+}
diff --git a/src/main.cpp b/src/main.cpp
index 75a842d..19ef8fd 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -11,6 +11,7 @@
  * */
 
 #include <cmath>
+#include <complex>
 #include <fstream>
 #include <functional>
 #include <omp.h>
@@ -29,6 +30,20 @@
 Particle p1(CHARGE, MASS, vec_3d{20., 0., 20.}, vec_3d{0., 25., 0.});
 Particle p2(CHARGE, MASS, vec_3d{25., 25., 0.}, vec_3d{0., 40., 5.});
 
+vec_3d analytical_solution_particle_1(double t)
+{
+    double w_0 = T / MASS;
+    double w_z2 = (50. * V / 1000.) / (MASS * 500. * 500.);
+    double w_p = (w_0 + std::sqrt(w_0 * w_0 - 2. * w_z2)) / 2.;
+    double w_n = (w_0 - std::sqrt(w_0 * w_0 - 2. * w_z2)) / 2.;
+    double A_p = (25. + w_n * 20.) / (w_n - w_p);
+    double A_n = -(25. + w_p * 20.) / (w_n - w_p);
+    std::complex<double> f = A_p * std::exp(std::complex<double>(0., -w_p * t)) +
+                             A_n * std::exp(std::complex<double>(0., -w_n * t));
+    vec_3d res{std::real(f), std::imag(f), 20. * std::cos(std::sqrt(w_z2) * t)};
+    return res;
+}
+
 void simulate_single_particle()
 {
     DEBUG("Inside single particle sim");
@@ -37,7 +52,8 @@ void simulate_single_particle()
     double time = 50.; // microseconds
 
     DEBUG("Write to dir");
-    trap.write_simulation_to_dir("output/simulate_single_particle", time, N, "rk4", false);
+    trap.write_simulation_to_dir("output/simulate_single_particle", time, N,
+                                 "rk4", false);
 }
 
 void simulate_two_particles()
@@ -56,22 +72,40 @@ void simulate_two_particles()
 void simulate_single_particle_with_different_steps()
 {
 
-    double time = 50; // microseconds
+    double time = 50.; // microseconds
+
+    std::ofstream ofile;
 
     for (int i = 0; i < 4; i++) {
-        int steps = 4000 * (i + 1);
+        int steps = 4000 * std::pow(2, i);
+        double dt = time / (double)steps;
+        std::string path = "output/relative_error/RK4/";
+        mkpath(path);
+        ofile.open(path + std::to_string(steps) + "_steps.txt");
         PenningTrap trap(std::vector<Particle>{p1});
-        trap.write_simulation_to_dir("output/N_steps/RK4/" +
-                                         std::to_string(steps) + "_steps",
-                                     time, steps, "rk4", false);
+        for (int i = 0; i < steps; i++) {
+            trap.evolve_RK4(dt);
+            ofile << arma::norm(trap.get_r(0) -
+                         analytical_solution_particle_1(trap.get_t()))
+                  << "\n";
+        }
+        ofile.close();
     }
 
     for (int i = 0; i < 4; i++) {
-        int steps = 4000 * (i + 1);
+        int steps = 4000 * std::pow(2, i);
+        double dt = time / (double)steps;
+        std::string path = "output/relative_error/euler/";
+        mkpath(path);
+        ofile.open(path + std::to_string(steps) + "_steps.txt");
         PenningTrap trap(std::vector<Particle>{p1});
-        trap.write_simulation_to_dir("output/N_steps/euler/" +
-                                         std::to_string(steps) + "_steps",
-                                     time, steps, "euler", false);
+        for (int i = 0; i < steps; i++) {
+            trap.evolve_forward_euler(dt);
+            ofile << arma::norm(trap.get_r(0) -
+                         analytical_solution_particle_1(trap.get_t()))
+                  << "\n";
+        }
+        ofile.close();
     }
 }
 
@@ -105,7 +139,7 @@ void simulate_100_particles_with_time_potential()
     std::ofstream ofile;
 
     double freq = freq_start;
-    for (size_t i=0; i<freq_iterations; i++) {
+    for (size_t i = 0; i < freq_iterations; i++) {
         res[0][i] = freq;
         freq += freq_increment;
     }
@@ -113,24 +147,23 @@ void simulate_100_particles_with_time_potential()
 #pragma omp parallel for collapse(2) num_threads(4)
     for (size_t i = 0; i < 3; i++) {
         for (size_t j = 0; j < freq_iterations; j++) {
-            PenningTrap trap((unsigned)100, T,
-                             std::bind([](double f, double r, double t) {
-                                 return (25. * V / 1000.) *
-                                        (1. +
-                                         f * std::cos(r * t));
-                             }, amplitudes[i], res[0][j], std::placeholders::_1),
-                             500., 0.);
+            PenningTrap trap(
+                (unsigned)100, T,
+                std::bind(
+                    [](double f, double r, double t) {
+                        return (25. * V / 1000.) * (1. + f * std::cos(r * t));
+                    },
+                    amplitudes[i], res[0][j], std::placeholders::_1),
+                500., 0.);
             res[i + 1][j] =
                 trap.fraction_of_particles_left(500., 40000, "rk4", false);
         }
     }
 
     ofile.open(path + "res.txt");
-    for (size_t i=0; i<freq_iterations; i++) {
-        ofile << res[0][i] << "," 
-            << res[1][i] << ","
-            << res[2][i] << ","
-            << res[3][i] << "\n";
+    for (size_t i = 0; i < freq_iterations; i++) {
+        ofile << res[0][i] << "," << res[1][i] << "," << res[2][i] << ","
+              << res[3][i] << "\n";
     }
     ofile.close();
 }
@@ -143,11 +176,11 @@ int main()
 
     simulate_two_particles();
 
-    //simulate_single_particle_with_different_steps();
+     simulate_single_particle_with_different_steps();
 
-    //simulate_100_particles();
+    // simulate_100_particles();
 
-    //simulate_100_particles_with_time_potential();
+    // simulate_100_particles_with_time_potential();
 
     double end = omp_get_wtime();
 
diff --git a/src/scripts/plot_2_particles.py b/src/scripts/plot_2_particles.py
new file mode 100644
index 0000000..17f327c
--- /dev/null
+++ b/src/scripts/plot_2_particles.py
@@ -0,0 +1,39 @@
+import matplotlib.pyplot as plt
+import numpy as np
+
+def main():
+    files = [
+        "output/simulate_2_particles/no_interaction/particle_0_r.txt",
+        "output/simulate_2_particles/no_interaction/particle_1_r.txt",
+        "output/simulate_2_particles/with_interaction/particle_0_r.txt",
+        "output/simulate_2_particles/with_interaction/particle_1_r.txt"
+    ]
+    labels = [
+        "particle 1 no interaction",
+        "particle 2 no interaction",
+        "particle 1 with interaction",
+        "particle 2 with interaction",
+    ]
+    colors = [
+        "lightskyblue",
+        "lightskyblue",
+        "salmon",
+        "salmon"
+    ]
+    for label, color, file in zip(labels, colors, files):
+        with open(file) as f:
+            lines = f.readlines()
+            t = np.linspace(0, 50, len(lines))
+            r = np.array([list(map(float, line.strip().split(","))) for line in lines])
+            plt.plot(r[:,0], r[:,1], label=label, color=color)
+    
+    plt.xlabel(r"x $(\mu m)$")
+    plt.ylabel(r"y $(\mu m)$")
+    plt.title(r"2 particles with and without interactions.")
+    # plt.legend()
+    # plt.show()
+    plt.savefig("../latex/images/plot_2_particles_xy.pdf")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/src/scripts/plot_3d.py b/src/scripts/plot_3d.py
new file mode 100644
index 0000000..811bf6a
--- /dev/null
+++ b/src/scripts/plot_3d.py
@@ -0,0 +1,40 @@
+import matplotlib.pyplot as plt
+import numpy as np
+
+def main():
+    files = [
+        "output/simulate_2_particles/no_interaction/particle_0_r.txt",
+        "output/simulate_2_particles/no_interaction/particle_1_r.txt",
+        "output/simulate_2_particles/with_interaction/particle_0_r.txt",
+        "output/simulate_2_particles/with_interaction/particle_1_r.txt"
+    ]
+    labels = [
+        "particle 1 no interaction",
+        "particle 2 no interaction",
+        "particle 1 with interaction",
+        "particle 2 with interaction",
+    ]
+    colors = [
+        "lightskyblue",
+        "deepskyblue",
+        "salmon",
+        "darkred"
+    ]
+    ax = plt.figure().add_subplot(projection="3d")
+    for label, color, file in zip(labels, colors, files):
+        with open(file) as f:
+            lines = f.readlines()
+            t = np.linspace(0, 50, len(lines))
+            r = np.array([list(map(float, line.strip().split(","))) for line in lines])
+            ax.plot(r[:,0], r[:,1], r[:,2], label=label, color=color)
+    
+    ax.set_xlabel(r"x $(\mu m)$")
+    ax.set_ylabel(r"y $(\mu m)$")
+    ax.set_zlabel(r"z $(\mu m)$")
+    plt.title(r"2 particles with and without interactions.")
+    plt.legend()
+    plt.savefig("../latex/images/3d_plot.pdf")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/src/scripts/plot_phase_space.py b/src/scripts/plot_phase_space.py
new file mode 100644
index 0000000..5a4c3ac
--- /dev/null
+++ b/src/scripts/plot_phase_space.py
@@ -0,0 +1,53 @@
+import matplotlib.pyplot as plt
+import numpy as np
+
+def main():
+    directories = {
+        "output/simulate_2_particles/no_interaction/",
+        "output/simulate_2_particles/with_interaction/",
+    }
+    titles = {
+        "particles without interaction",
+        "particles with interaction"
+    }
+    files = [
+        "particle_0_r.txt",
+        "particle_0_v.txt",
+        "particle_1_r.txt",
+        "particle_1_v.txt",
+    ]
+    labels = [
+        r"particle 1 r",
+        r"particle 1 v",
+        r"particle 2 r",
+        r"particle 2 v",
+    ]
+    colors = [
+        "lightskyblue",
+        "deepskyblue",
+        "salmon",
+        "tomato",
+    ]
+    fig1, axs1 = plt.subplots(2,1)
+    fig2, axs2 = plt.subplots(2,1)
+    for i, (dir, title) in enumerate(zip(directories, titles)):
+        for label, color, file in zip(labels, colors, files):
+            with open(dir+file) as f:
+                lines = f.readlines()
+                t = np.linspace(0, 50, len(lines))
+                r = np.array([list(map(float, line.strip().split(","))) for line in lines])
+                axs1[i].plot(t, r[:,0], label=label, color=color)
+                axs2[i].plot(t, r[:,2], label=label, color=color)
+        
+            axs1[i].set(xlabel=r"t $(\mu s)$", ylabel = r"z $(\mu m)$")
+        
+        axs1[i].legend()
+        axs1[i].set_title(title)
+
+    # plt.show()
+    fig1.savefig("../latex/images/phase_space_2_particles_x.pdf")
+    fig2.savefig("../latex/images/phase_space_2_particles_z.pdf")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/src/scripts/plot_relative_error.py b/src/scripts/plot_relative_error.py
new file mode 100644
index 0000000..54214c8
--- /dev/null
+++ b/src/scripts/plot_relative_error.py
@@ -0,0 +1,50 @@
+import matplotlib.pyplot as plt
+import numpy as np
+
+def main():
+    directories = {
+        "output/relative_error/RK4/",
+        "output/relative_error/euler/",
+    }
+    files = [
+        "4000_steps.txt",
+        "8000_steps.txt",
+        "16000_steps.txt",
+        "32000_steps.txt",
+    ]
+    labels = [
+        r"4000 steps",
+        r"8000 steps",
+        r"16000 steps",
+        r"32000 steps",
+    ]
+    titles = [
+        "Relative error for the RK4 method",
+        "Relative error for the forward Euler method"
+    ]
+    fig1, axs1 = plt.subplots(2,1)
+    for i, (dir, title) in enumerate(zip(directories, titles)):
+        max_err = []
+        for label, file in zip(labels, files):
+            with open(dir+file) as f:
+                lines = f.readlines()
+                t = np.linspace(0, 50, len(lines))
+                r = np.array([float(line.strip()) for line in lines])
+                max_err.append(max(r))
+                axs1[i].plot(t, r, label=label)
+        
+            axs1[i].set(xlabel=r"t $(\mu s)$", ylabel = r"relative_error $(\mu m)$")
+
+        
+        axs1[i].legend()
+        axs1[i].set_title(title)
+
+        conv_rate = 1/3 * sum([np.log10(max_err[i+1]/max_err[i])/np.log10(.5) for i in range(3)])
+        print(conv_rate)
+
+    plt.show()
+    # fig1.savefig("../latex/images/phase_space_2_particles_x.pdf")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/src/scripts/plot_single_particle.py b/src/scripts/plot_single_particle.py
index 2874030..f4bdc6a 100644
--- a/src/scripts/plot_single_particle.py
+++ b/src/scripts/plot_single_particle.py
@@ -22,7 +22,8 @@ def main():
     plt.ylabel(r"z $(\mu m)$")
     plt.title(r"Movement of a single particle in the x direction")
     plt.legend()
-    plt.savefig("../latex/images/single_particle.pdf")
+    # plt.savefig("../latex/images/single_particle.pdf")
+    plt.show()
 
 
 if __name__ == "__main__":