diff --git a/latex/main.pdf b/latex/main.pdf
index 6e3a093..979dd80 100644
Binary files a/latex/main.pdf and b/latex/main.pdf differ
diff --git a/latex/references/references.bib b/latex/references/references.bib
index 236bbd0..38b1079 100644
--- a/latex/references/references.bib
+++ b/latex/references/references.bib
@@ -7,3 +7,12 @@
 	publisher= "CRC Press",
 	year	= "2016",
 }
+
+@article{rk4_method,
+    author = "Morten Hjorth-Jensen",
+    title = "Computational Physics, Lecture Notes Fall 2015",
+    journal = "Department of Physics, University of Oslo",
+    year = "2015",
+    chapter = "8.4",
+    pages = "250--252",
+}
diff --git a/latex/sections/methods.tex b/latex/sections/methods.tex
index fb2b32c..80ca419 100644
--- a/latex/sections/methods.tex
+++ b/latex/sections/methods.tex
@@ -6,18 +6,31 @@
 
 \section{Methods}
 
-Newton's second law of motion is defined as
+For a multi-particles system, we need to modify $\vb{F}$ to account for the 
+force of other particles in the system acting upon each other. To do that, we
+add another term to $\vb{F}$
 \begin{equation}
-    \frac{d^2\vb{r}}{dt^2} = \frac{\vb{F}\left(t, \frac{d\vb{r}}{dt}, \vb{r}\right)}{m}.
+    \vb{F}_i(t, \vb{v}_i, \vb{r}_i) = q_i \vb{E}(t, \vb{r}_i) + q_i \vb{v}_i \cross \vb{B} - \vb{E}_p(t, \vb{r}_i),
+\end{equation}
+where $i$ and $j$ are particle indices and
+\begin{equation}
+    \vb{E}_p(t, \vb{r}_i) = q_i k_e \sum_{j \neq i} 
+    q_j \frac{\vb{r_i} - \vb{r_j}}{\left| \vb{r_i} - \vb{r_j} \right|^3}.
+    \label{eq:}
+\end{equation}
+
+Newton's second law for a particle $i$ is then
+\begin{equation}
+    \frac{d^2\vb{r}_i}{dt^2} = \frac{\vb{F}\left(t, \frac{d\vb{r}_i}{dt}, \vb{r_i}\right)}{m_i},
     \label{eq:newtonlaw2}
 \end{equation}
 We can then rewrite the second order ODE from equation~\ref{eq:newtonlaw2} 
 into a set of coupled first order ODEs. 
-We now redefine Newton's second law of motion as
+We now rewrite Newton's second law of motion as
 \begin{equation}
     \begin{split}
-        \frac{d\vb{r}}{dt} &= \vb{v} \\
-        \frac{d\vb{v}}{dt} &= \frac{\vb{F}(t, \vb{v}, \vb{r})}{m}.
+        \frac{d\vb{r}_i}{dt} &= \vb{v}_i \\
+        \frac{d\vb{v}_i}{dt} &= \frac{\vb{F}(t, \vb{v}_i, \vb{r}_i)}{m_i}.
     \end{split}
     \label{eq:coupled}
 \end{equation}
@@ -33,9 +46,8 @@ expressed as
         \vb{v}_{i+1} &= \vb{v}_i + h \cdot \frac{\vb{v}_i}{dt} = \vb{v}_i + h \cdot \frac{\vb{F}\left( t_i, \vb{v}_i, \vb{r}_i \right)}{m},
     \end{split}
 \end{equation}
-where $\vb{v}_i$ is the current velocity of the particle, $\vb{r}_i$ is the
-current position of the particle, $m$ is the mass of the particle, and $h$
-is a predetermined timestep.
+where $i$ is the current time step of the particle $m$ is the mass of the particle, and $h$
+is the time step length.
 
 When dealing with a multi-particle system, we need to ensure that we do not
 update the position of any particles until every particle has calculated their
@@ -50,12 +62,12 @@ Algorithm~\ref{algo:forwardeuler} provides an overview on how that can be achiev
     \label{algo:forwardeuler}
         \begin{algorithmic}
             \Procedure{Evolve forward Euler}{$particles, dt$}
-            \State $new\_state \leftarrow particles$ \Comment{Create a copy of the particles}
-            \For{ $i = 1, 2, \ldots , \left| particles \right|$ }
-            \State $new\_state_i.\vb{v} \leftarrow particles_i.\vb{v} + dt \cdot \frac{\vb{F}}{m}$
-            \State $new\_state_i.\vb{r} \leftarrow particles_i.\vb{r} + dt \cdot particles_i.\vb{v}$
+            \State $N \leftarrow \text{Number of particles in } particles$
+            \State $a \leftarrow \text{Calculate } \frac{\vb{F_i}}{m_i} \text{ for each particle in } particles$
+            \For{ $i = 1, 2, \ldots , N$ }
+            \State $particles_i.\vb{r} \leftarrow particles_i.\vb{r} + dt \cdot particles_i.\vb{v}$
+            \State $particles_i.\vb{v} \leftarrow particles_i.\vb{v} + dt \cdot a_i$
             \EndFor
-            \State $particles \leftarrow new\_state$
             \EndProcedure
         \end{algorithmic}
     \end{algorithm}
@@ -88,6 +100,19 @@ where
     \end{split}
 \end{equation}
 
+In order to find each $\vb{k}_{\vb{r},i}$ and $\vb{k}_{\vb{v},i}$, 
+we need to first compute all $\vb{k}_{\vb{r},i}$ and $\vb{k}_{\vb{v},i}$ 
+for all particles, then we can update the particles in order to compute 
+$\vb{k}_{\vb{r},i+1}$ and $\vb{k}_{\vb{v},i+1}$. In order for the algorithm
+to work, we need to save a copy of each particle before starting so that we
+can update the particles correctly for each step.
+
+This approach would require 8 loops to be able to complete the calculation since
+we cannot update the particles until after all $\vb{k}$ values have been 
+computed, however if we create a temporary array that holds particles, we can
+put the updated particles in there, and then use that array in the next loop.
+around 8 separate loops to calculate all $\vb{k}$ values
+
 \newpage
 
 \begin{figure}