diff --git a/charge/charge.c b/charge/charge.c
index 2158dfdc2b20123b2103dcc7b6e462ba151edaf3..9bc963e2edf402da71122541c57413d429934eb0 100644
--- a/charge/charge.c
+++ b/charge/charge.c
@@ -75,12 +75,13 @@ double random_charge_q(bool positive) {
     return (rand() % 100) * (positive ? 1 : -1);
 }
 
+int sign(double a) { return (a < 0) ? -1 : 1; }
+
 vec2 electristatic_force(charge_t a, charge_t b) {
+    bool attractive = sign(a.q) != sign(b.q);
     vec2 rab = vec2_sub(b.pos, a.pos);
     vec2 res = vec2_mul(
-        rab, K * (((a.q * 10e-6) * (b.q * 10e-6)) / pow(vec2_norm(rab), 2)));
-
-    printf("COEF: %.25f\n", (a.q * 10e-6) * (b.q * 10e-6));
+        rab, -K * (((a.q * 10e-6) * (b.q * 10e-6)) / pow(vec2_norm(rab), 2)));
     return res;
 }
 
@@ -100,18 +101,36 @@ vec2 resulting_electrostatic_force(charge_t *charges, int num_charges,
     return resulting;
 }
 
-void update_charge_pos(charge_t *charges, int num_charges, charge_t *a) {
-    vec2 resulting = resulting_electrostatic_force(charges, num_charges, *a);
-    a->pos = vec2_add(a->pos, resulting);
+vec2 compute_initial_vecocity(charge_t c, charge_t *charges, int num_charges) {
+    double norm = 0;
+    double angle = rand() % 360;
+    return vec2_create(norm * cos(angle), norm * sin(angle));
+}
+
+vec2 compute_acceleration(charge_t c, charge_t *charges, int num_charges) {
+    return vec2_mul(resulting_electrostatic_force(charges, num_charges, c),
+                    1.0 / CHARGE_MASS);
+}
+
+vec2 compute_first_pos(charge_t c, charge_t *charges, int num_charges) {
+    vec2 velocity = compute_initial_vecocity(c, charges, num_charges);
+    vec2 a = compute_acceleration(c, charges, num_charges);
+    vec2 pos = vec2_add(vec2_add(c.pos, vec2_mul(velocity, DELTA_T)),
+                        vec2_mul(a, pow(DELTA_T, 2) / 2));
+    return pos;
+}
+
+vec2 compute_next_pos(charge_t c, charge_t *charges, int num_charges) {
+    vec2 a = compute_acceleration(c, charges, num_charges);
+    vec2 pos = vec2_add(vec2_sub(vec2_mul(c.pos, 2), c.prev_pos),
+                        vec2_mul(a, pow(DELTA_T, 2)));
+    return pos;
 }
 
 void update_simulation(charge_t *charges, int num_charges) {
     for (int i = 0; i < num_charges; i++) {
-        vec2 resulting =
-            resulting_electrostatic_force(charges, num_charges, charges[i]);
-        printf("RES: %.19f\n", vec2_norm(resulting));
-        charges[i].pos = vec2_add(charges[i].pos, resulting);
+        vec2 tmp = charges[i].pos;
+        charges[i].pos = compute_next_pos(charges[i], charges, num_charges);
+        charges[i].prev_pos = tmp;
     }
-
-    printf("%.19f, %.19f\n", charges[0].pos.x, charges[0].pos.x);
 }
\ No newline at end of file
diff --git a/charge/charge.h b/charge/charge.h
index 57fc6687f5111c5b5e6d10e1001fcdf9a975666a..8a3c229fe0a584b32f81b7e89b1382d6e3942ddf 100644
--- a/charge/charge.h
+++ b/charge/charge.h
@@ -23,6 +23,8 @@
 #define ELEMENTARY_CHARGE 1.602176634e-19
 #define CHARGE_RADIUS 10
 #define EPS 0.015
+#define DELTA_T 10e-4
+#define CHARGE_MASS 1500
 #define UNIV_MIN_X 0
 #define UNIV_MAX_X 1
 #define UNIV_MIN_Y 0
@@ -31,6 +33,7 @@
 typedef struct _charge_t {
     double q;
     vec2 pos;
+    vec2 prev_pos;
 } charge_t;
 
 /**
@@ -140,14 +143,14 @@ vec2 resulting_electrostatic_force(charge_t *charges, int num_charges,
                                    charge_t a);
 
 /**
- * @brief Update the pos of the given charge based on the resulting
- * electrostatic force
+ * @brief Compute the first pos for the dynamic simulation
  *
+ * @param c
  * @param charges
  * @param num_charges
- * @param a
+ * @return vec2
  */
-void update_charge_pos(charge_t *charges, int num_charges, charge_t *a);
+vec2 compute_first_pos(charge_t c, charge_t *charges, int num_charges);
 
 /**
  * @brief Update the whole simulation
diff --git a/draw/draw.c b/draw/draw.c
index 021825507bdf84498ec63dcfde193fdc065f6ab8..5ac085d8de3ecbb72bb22b7d325f6283ea9bfce0 100644
--- a/draw/draw.c
+++ b/draw/draw.c
@@ -49,12 +49,19 @@ void draw_line(vec2 start, vec2 end) {
 }
 
 void get_heatmap_color(double value, double *red, double *green, double *blue) {
-    double aR = 0.0;
-    double aG = 0.0;
-    double aB = 0.0; // RGB for our 1st color (blue in this case).
-    double bR = 0.31;
-    double bG = 0.23;
-    double bB = 0.31; // RGB for our 2nd color (red in this case).
+    double aR = 1.0;
+    double aG = 1.0;
+    double aB = 1.0; // RGB for our 1st color
+    double bR = 0.67;
+    double bG = 0.12;
+    double bB = 0.47; // RGB for our 2nd color
+
+    // if (value < 0.3) {
+    //     *red = 1.0;
+    //     *green = 1.0;
+    //     *blue = 1.0;
+    //     return;
+    // }
 
     *red = (bR - aR) * value + aR;   // Evaluated as -255*value + 255.
     *green = (bG - aG) * value + aG; // Evaluates as 0.
@@ -170,12 +177,14 @@ void draw_vector_field(charge_t *charges, int num_charges, vec2 *points,
     for (int i = 0; i < nb_points; i++) {
         vec2 a = position_to_coordinates(w, h, points[i]);
 
+        glLineWidth(2);
         glBegin(GL_LINES);
-        glColor3f(0.61, 0.87, 0.61);
+        glColor3f(0.18, 0.46, 0.18);
         glVertex2d(a.x, a.y);
-        a = vec2_add(a, vec2_mul(vec2_div(vs[i], max), 10));
+        a = vec2_add(a, vec2_mul(vec2_div(vs[i], max), 40));
         glVertex2d(a.x, a.y);
         glEnd();
+        glLineWidth(1);
     }
 }
 
diff --git a/main.c b/main.c
index 4ff31394bf491ddcd7cee718ce4001cd759db475..1b37a23db44c29ec1066862dba5043c563fb26e4 100644
--- a/main.c
+++ b/main.c
@@ -30,12 +30,12 @@
 #endif
 
 #define ONE_SECOND_IN_MILLISECONDS 1000
-#define REFRESH_RATE 200
+#define REFRESH_RATE 60
 
 #define WINDOW_WIDTH 800
 #define WINDOW_HEIGHT 800
 #define NB_CHARGES 5
-#define GRID_SPACING 10
+#define GRID_SPACING 25
 #define NB_POINTS                                                              \
     ((WINDOW_WIDTH / GRID_SPACING) * (WINDOW_WIDTH / GRID_SPACING))
 #define MAXIMUM_CHARGES 30 // Can be replace by a dynamic array
@@ -47,7 +47,7 @@ vec2 default_charges_pos[MAXIMUM_CHARGES];
 double default_charges[MAXIMUM_CHARGES];
 charge_t *charges;
 
-bool toggle_lines_view = false;
+bool toggle_lines_view = true;
 bool toggle_vector_view = false;
 bool toggle_heatmap_view = false;
 
@@ -85,7 +85,7 @@ void generate_points(vec2 points[]) {
  *
  */
 void display() {
-    glClearColor(0, 0, 0, 1);
+    glClearColor(1, 1, 1, 1);
     glClear(GL_COLOR_BUFFER_BIT);
 
     if (toggle_lines_view) {
@@ -186,8 +186,11 @@ int main(int argc, char **argv) {
     }
 
     // Init charges
-    for (int i = 0; i < charge_count; i++)
+    for (int i = 0; i < charge_count; i++) {
         charges[i] = charge_create(default_charges[i], default_charges_pos[i]);
+        charges[i].prev_pos = charges[i].pos;
+        charges[i].pos = compute_first_pos(charges[i], charges, charge_count);
+    }
 
     // Generate all field lines starting points
     generate_points(starting_points);
@@ -209,9 +212,8 @@ int main(int argc, char **argv) {
     glutMouseFunc(handle_mouse_input);
     glutKeyboardFunc(keyboard);
 
-    // glutTimerFunc(ONE_SECOND_IN_MILLISECONDS / REFRESH_RATE, update_timer,
-    // 0); glutTimerFunc(ONE_SECOND_IN_MILLISECONDS / REFRESH_RATE, draw_timer,
-    // 0);
+    glutTimerFunc(ONE_SECOND_IN_MILLISECONDS / REFRESH_RATE, update_timer, 0);
+    glutTimerFunc(ONE_SECOND_IN_MILLISECONDS / REFRESH_RATE, draw_timer, 0);
 
     glutMainLoop();