diff --git a/nbody_problem/cpp/main.cpp b/nbody_problem/cpp/main.cpp
new file mode 100644
index 0000000..ed51c16
--- /dev/null
+++ b/nbody_problem/cpp/main.cpp
@@ -0,0 +1,170 @@
+#include <iostream>
+#include <math.h>
+#include <numeric>
+#include <vector>
+
+typedef double real;
+using namespace std;
+
+class Star {
+public:
+  real m;
+  vector<real> r;
+  vector<real> v;
+  vector<real> a, a0;
+
+  Star() {
+    r.assign(3, 0);
+    v.assign(3, 0);
+    a.assign(3, 0);
+    a0.assign(3, 0);
+  }
+
+  Star(real mass, vector<real> pos, vector<real> vel) {
+    m = mass;
+    r = pos;
+    v = vel;
+  }
+
+  friend ostream &operator<<(ostream &so, const Star &si) {
+    so << si.m << " " << si.r[0] << " " << si.r[1] << " " << si.r[2] << " "
+       << si.v[0] << " " << si.v[1] << " " << si.v[2] << endl;
+    return so;
+  }
+};
+
+class Cluster : public Star {
+protected:
+public:
+  vector<Star> s;
+  Cluster() : Star() {}
+
+  void acceleration() {
+    for (vector<Star>::iterator si = s.begin(); si != s.end(); ++si)
+      si->a.assign(3, 0);
+
+    for (vector<Star>::iterator si = s.begin(); si != s.end(); ++si) {
+      vector<real> rij(3);
+      real init = 0.0;
+
+      for (vector<Star>::iterator sj = s.begin(); sj != s.end(); ++sj) {
+        if (si != sj) {
+
+          for (int i = 0; i != 3; ++i)
+            rij[i] = si->r[i] - sj->r[i];
+
+          real RdotR = inner_product(rij.begin(), rij.end(), rij.begin(), init);
+          real apre = 1. / sqrt(RdotR * RdotR * RdotR);
+
+          for (int i = 0; i != 3; ++i) {
+            si->a[i] = sj->m * -1 * apre * rij[i];
+          }
+        }
+      }
+    }
+  }
+
+  void updatePositions(real dt) {
+    for (vector<Star>::iterator si = s.begin(); si != s.end(); ++si) {
+      si->a0 = si->a;
+      for (int i = 0; i != 3; ++i)
+        si->r[i] += dt * si->v[i] + 0.5 * dt * dt * si->a0[i];
+    }
+  }
+
+  void updateVelocities(real dt) {
+
+    for (vector<Star>::iterator si = s.begin(); si != s.end(); ++si) {
+      for (int i = 0; i != 3; ++i)
+        si->v[i] += 0.5 * dt * (si->a0[i] + si->a[i]);
+      si->a0 = si->a;
+    }
+  }
+
+  vector<real> energies() {
+    real init = 0;
+    vector<real> E(3), rij(3);
+    E.assign(3, 0);
+
+    for (vector<Star>::iterator si = s.begin(); si != s.end(); ++si)
+      E[1] += 0.5 * si->m *
+              inner_product(si->v.begin(), si->v.end(), si->v.begin(), init);
+
+    // Potential energy
+    for (vector<Star>::iterator si = s.begin(); si != s.end(); ++si) {
+      for (vector<Star>::iterator sj = si + 1; sj != s.end(); ++sj) {
+        for (int i = 0; i != 3; ++i)
+          rij[i] = si->r[i] - sj->r[i];
+        E[2] -= si->m * sj->m /
+                sqrt(inner_product(rij.begin(), rij.end(), rij.begin(), init));
+      }
+    }
+    E[0] = E[1] + E[2];
+    return E;
+  }
+
+  // Print function
+  friend ostream &operator<<(ostream &so, Cluster &cl) {
+    for (vector<Star>::iterator si = cl.s.begin(); si != cl.s.end(); ++si)
+      so << *si;
+    return so;
+  }
+};
+
+int main(int argc, char* argv[]) {
+
+  Cluster cl;
+  real m;
+  int dummy;
+  vector<real> r(3), v(3);
+
+  // Read input data from the command line (makeplummer | dumbp)
+  do {
+    cin >> dummy;
+    cin >> m;
+    for (int i = 0; i != 3; ++i)
+      cin >> r[i];
+    for (int i = 0; i != 3; ++i)
+      cin >> v[i];
+    cl.s.push_back(Star(m, r, v));
+  } while (!cin.eof());
+
+  cl.s.pop_back();
+
+  // Compute initial energu of the system
+  vector<real> E(3), E0(3);
+  E0 = cl.energies();
+
+  // Start time, end time and simulation step
+  real t = 0.0;
+  real tend;
+
+  if (argc > 1)
+    tend = strtod(argv[1], NULL);
+  else
+    tend = 10.0;
+
+  real dt = 1e-3;
+  int k = 0;
+
+  for (int i = 0; i < 50; i++) {
+    cl.acceleration();
+
+    while (t < tend) {
+      cl.updatePositions(dt);
+
+      cl.acceleration();
+
+      cl.updateVelocities(dt);
+
+      t += dt;
+      k += 1;
+      if (k % 100 == 0) {
+        E = cl.energies();
+      }
+    }
+    t = 0;
+  }
+
+  return 0;
+}
diff --git a/nbody_problem/python/compiled_methods.py b/nbody_problem/python/compiled_methods.py
new file mode 100644
index 0000000..7cad8ab
--- /dev/null
+++ b/nbody_problem/python/compiled_methods.py
@@ -0,0 +1,108 @@
+import sys
+import math
+import timeit
+
+import numpy as np
+import pandas as pd
+
+from transonic import jit
+
+def load_input_data(path):
+
+    df = pd.read_csv(
+        path, names = ["mass", "x", "y", "z", "vx", "vy", "vz"], delimiter=r"\s+"
+    )
+
+    masses = df["mass"].values.copy()
+    positions = df.loc[:, ["x", "y", "z"]].values.copy()
+    velocities = df.loc[:, ["vx", "vy", "vz"]].values.copy()
+
+    return masses, positions, velocities
+
+def compute_accelerations(accelerations, masses, positions):
+    nb_particles = masses.size
+
+    for index_p0 in range(nb_particles - 1):
+        position0 = positions[index_p0]
+        mass0 = masses[index_p0]
+        
+        # TODO: Use compiled methods like Numba & Pythran in vectorized approach.
+        # Issue: https://github.com/khushi-411/numpy-benchmarks/issues/4
+        for index_p1 in range(index_p0 + 1, nb_particles):
+            mass1 = masses[index_p1]
+
+            vectors = position0 - positions[index_p1]
+
+            distances = (vectors**2).sum()
+            coefs = 1./distances**1.5
+
+            accelerations[index_p0] += mass1 * -1 * vectors * coefs
+            accelerations[index_p1] += mass0 * vectors * coefs
+
+    return accelerations
+
+@jit
+def loop(
+        time_step, nb_steps, masses, positions,
+        velocities):
+
+    accelerations = np.zeros_like(positions)
+    accelerations1 = np.zeros_like(positions)
+
+    accelerations = compute_accelerations(accelerations, masses, positions)
+
+    time = 0.0
+    energy0, _, _ = compute_energies(masses, positions, velocities)
+    energy_previous = energy0
+
+    for step in range(nb_steps):
+        positions += time_step*velocities + 0.5*accelerations*time_step**2
+
+        accelerations, accelerations1 = accelerations1, accelerations
+        accelerations.fill(0)
+        accelerations = compute_accelerations(accelerations, masses, positions)
+
+        velocities += 0.5*time_step*(accelerations+accelerations1)
+
+        time += time_step
+
+        if not step % 100:
+            energy, _, _ = compute_energies(masses, positions, velocities)
+            energy_previous = energy
+
+    return energy, energy0
+
+def compute_energies(masses, positions, velocities):
+
+    ke = 0.5 * masses @ np.sum(velocities**2, axis=1)
+
+    nb_particles = masses.size
+    pe = 0.0
+    for index_p0 in range(nb_particles - 1):
+
+        mass0 = masses[index_p0]
+        for index_p1 in range(index_p0 + 1, nb_particles):
+
+            mass1 = masses[index_p1]
+            vector = positions[index_p0] - positions[index_p1]
+
+            distance = math.sqrt((vector**2).sum())
+
+            pe -= (mass0*mass1) / distance**2
+
+    return ke+pe, ke, pe
+
+if __name__ == "__main__":
+
+    try:
+        time_end = float(sys.argv[2])
+    except IndexError:
+        time_end = 10.0
+
+    time_step = 0.001
+    nb_steps = int(time_end/time_step) + 1
+
+    path_input = sys.argv[1]
+    masses, positions, velocities = load_input_data(path_input)
+
+    print('time taken:', timeit.timeit('loop(time_step, nb_steps, masses, positions, velocities)', globals=globals(), number=50))
diff --git a/nbody_problem/python/numpy_python.py b/nbody_problem/python/numpy_python.py
new file mode 100644
index 0000000..41bd461
--- /dev/null
+++ b/nbody_problem/python/numpy_python.py
@@ -0,0 +1,105 @@
+import sys
+import math
+import timeit
+
+import numpy as np
+import pandas as pd
+
+def load_input_data(path):
+
+    df = pd.read_csv(
+        path, names = ["mass", "x", "y", "z", "vx", "vy", "vz"], delimiter=r"\s+"
+    )
+
+    masses = df["mass"].values.copy()
+    positions = df.loc[:, ["x", "y", "z"]].values.copy()
+    velocities = df.loc[:, ["vx", "vy", "vz"]].values.copy()
+
+    return masses, positions, velocities
+
+def compute_accelerations(accelerations, masses, positions):
+    nb_particles = masses.size
+
+    for index_p0 in range(nb_particles - 1):
+        position0 = positions[index_p0]
+        mass0 = masses[index_p0]
+
+        vectors = position0 - positions[index_p0 + 1: nb_particles]
+
+        distances = (vectors**2).sum(axis=1)
+        coefs = 1./distances**1.5
+
+        accelerations[index_p0] += np.sum((masses[index_p0 + 1: nb_particles] * -1 * vectors.T * coefs).T, axis=0)
+        accelerations[index_p0 + 1: nb_particles] += (mass0 * vectors.T * coefs).T
+
+    return accelerations
+
+def numpy_loop(
+        time_step: float,
+        nb_steps: int,
+        masses: "float[]",
+        positions: "float[:,:]",
+        velocities: "float[:,:]",
+    ):
+
+    accelerations = np.zeros_like(positions)
+    accelerations1 = np.zeros_like(positions)
+
+    accelerations = compute_accelerations(accelerations, masses, positions)
+
+    time = 0.0
+
+    energy0, _, _ = compute_energies(masses, positions, velocities)
+    energy_previous = energy0
+
+    for step in range(nb_steps):
+        positions += time_step*velocities + 0.5*accelerations*time_step**2
+
+        accelerations, accelerations1 = accelerations1, accelerations
+        accelerations.fill(0)
+        accelerations = compute_accelerations(accelerations, masses, positions)
+
+        velocities += 0.5*time_step*(accelerations+accelerations1)
+
+        time += time_step
+
+        if not step%100:
+            energy, _, _ = compute_energies(masses, positions, velocities)
+            energy_previous = energy
+
+    return energy, energy0
+
+def compute_energies(masses, positions, velocities):
+
+    ke = 0.5 * masses @ np.sum(velocities**2, axis=1)
+
+    nb_particles = masses.size
+    pe = 0.0
+    for index_p0 in range(nb_particles - 1):
+
+        mass0 = masses[index_p0]
+        for index_p1 in range(index_p0 + 1, nb_particles):
+
+            mass1 = masses[index_p1]
+            vector = positions[index_p0] - positions[index_p1]
+
+            distance = math.sqrt((vector**2).sum())
+
+            pe -= (mass0*mass1) / distance**2
+
+    return ke+pe, ke, pe
+
+if __name__ == "__main__":
+
+    try:
+        time_end = float(sys.argv[2])
+    except IndexError:
+        time_end = 10.0
+
+    time_step = 0.001
+    nb_steps = int(time_end/time_step) + 1
+
+    path_input = sys.argv[1]
+    masses, positions, velocities = load_input_data(path_input)
+
+    print('time taken:', timeit.timeit('numpy_loop(time_step, nb_steps, masses, positions, velocities)', globals=globals(), number=50))
diff --git a/nbody_problem/python/plot.py b/nbody_problem/python/plot.py
new file mode 100644
index 0000000..7f67de3
--- /dev/null
+++ b/nbody_problem/python/plot.py
@@ -0,0 +1,74 @@
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+
+def plot(x, labels, list_df, names):
+
+    plt.figure(figsize=(350, 220), dpi=4)
+    plt.subplot(2, 1, 1)
+
+    width = 0.25
+    rect = [0, 0, 0]
+    colors = ['#4D77CF', 'navy']
+
+    list1 = []
+    for ind, list in enumerate(list_df):
+        if ind < 2:
+            list1.append(list)
+
+    for ind, time_taken in enumerate(list1):
+        rect[ind] = plt.bar(x + width*ind, time_taken, width, align='center', color=colors[ind % len(colors)], label=labels[ind])
+        plt.bar_label(rect[ind], padding=3, fontsize=250)
+ 
+    plt.xticks(x + width/2, names, fontsize=300)
+    plt.tick_params(axis='both', which='major', pad=60)
+    plt.yticks(fontsize=300)
+    plt.legend(fontsize=250)
+    
+    plt.ylabel(r'$\frac{Time\ (in\ secs)}{nParticles^{2}}$', fontsize=350)
+    plt.xlabel(r"$Number\ of\ Particles\ Simulated(nParticles)$", fontsize=270)
+    plt.title(r"$Library\ based\ Implementation$", fontsize=300)
+
+    plt.subplot(2, 1, 2)
+    rect = [0, 0, 0]
+    colors = ['#044F88', '#4DABCF', '#013243']
+
+    list2 = []
+    for ind, list in enumerate(list_df):
+        if ind >= 2:
+            list2.append(list)
+
+    j = 2
+    for ind, time_taken in enumerate(list2):
+        rect[ind] = plt.bar(x + width*ind, time_taken, width, align='center', color=colors[ind % len(colors)], label=labels[j])
+        j += 1
+        plt.bar_label(rect[ind], padding=3, fontsize=250)
+ 
+    plt.xticks(x + width, names, fontsize=300)
+    plt.yticks(fontsize=300)
+    plt.legend(fontsize=250)
+
+    plt.ylabel(r'$\frac{Time\ (in\ secs)}{nParticles^{2}}$', fontsize=350)
+    plt.tick_params(axis='both', which='major', pad=60)
+    plt.xlabel(r"$Number\ of\ Particles\ Simulated(nParticles)$", fontsize=270)
+    plt.title(r"$Compiler\ based\ Implementation$", fontsize=300)
+    
+    plt.tight_layout(pad=70.0)
+    plt.savefig("performance_benchmarking")
+
+if __name__ == "__main__":
+
+    data_path = "benchmarks/data/table_new_norm.csv"
+    df = pd.read_csv(data_path)
+    df = df.drop(['Unnamed: 0'], axis=1)
+    df = df.T
+    df.columns = ["NumPy", "Python", "C++", "Numba", "Pythran"]
+
+    labels = df.columns
+    names = ["32", "64", "128", "256"]
+
+    df = df.T
+    x = np.arange(len(df.columns))
+    list_df = df[0:].to_numpy().tolist()
+
+    plot(x, labels, list_df, names)
diff --git a/nbody_problem/python/pure_python.py b/nbody_problem/python/pure_python.py
new file mode 100644
index 0000000..0729e2c
--- /dev/null
+++ b/nbody_problem/python/pure_python.py
@@ -0,0 +1,154 @@
+import sys
+import math
+import timeit
+
+import numpy as np
+import pandas as pd
+
+def load_input_data(path):
+
+    df = pd.read_csv(
+        path, names = ["mass", "x", "y", "z", "vx", "vy", "vz"], delimiter=r"\s+"
+    )
+
+    masses = df["mass"].values.copy()
+    positions = df.loc[:, ["x", "y", "z"]].values.copy()
+    velocities = df.loc[:, ["vx", "vy", "vz"]].values.copy()
+
+    return masses, positions, velocities
+
+def compute_accelerations(accelerations, masses, positions):
+    nb_particles = len(masses)
+
+    for index_p0 in range(nb_particles - 1):
+        position0 = positions[index_p0]
+        mass0 = masses[index_p0]
+
+        for index_p1 in range(index_p0 + 1, nb_particles):
+            mass1 = masses[index_p1]
+
+            vector = [p0 - p1 for (p0, p1) in zip(position0, positions[index_p1])]
+
+            distance = 0
+            for i in vector:
+                distance += i**2
+            coefs = 1./distance**1.5
+
+            accelerations[index_p0] = [sum(i) for i in zip([vec_val * mass1 * -1 * coefs for vec_val in vector], accelerations[index_p0])]
+            accelerations[index_p1] = [sum(i) for i in zip([vec_val * mass0 * coefs for vec_val in vector], accelerations[index_p1])]
+
+    return accelerations
+
+def loop(
+    time_step: float,
+    nb_steps: int,
+    masses: "float[]",
+    positions: "float[:,:]",
+    velocities: "float[:,:]",
+):
+
+    accelerations = [[0.0 for _ in range(3)] for _ in range(len(positions))]
+    accelerations1 = [[0.0 for _ in range(3)] for _ in range(len(positions))]
+
+    accelerations = compute_accelerations(accelerations, masses, positions)
+
+    time = 0.0
+    energy0, _, _ = compute_energies(masses, positions, velocities)
+    energy_previous = energy0
+
+    for step in range(nb_steps):
+        pos_val = []
+        for (vel, acc) in zip(velocities, accelerations):
+            pos = []
+            for (v, a) in zip(vel, acc):
+                pos.append(v*time_step + 0.5*a*time_step**2)
+            pos_val.append(pos)
+
+        positions1 = []
+        for (pos_val1, pos1) in zip(pos_val, positions):
+            pos_1 = []
+            for (pos_val0, pos0) in zip(pos_val1, pos1):
+                pos_1.append(pos_val0 + pos0)
+            positions1.append(pos_1)
+        positions = positions1
+
+        accelerations, accelerations1 = accelerations1, accelerations
+        accelerations = [[0.0 for acc0 in acc1] for acc1 in accelerations]
+        accelerations = compute_accelerations(accelerations, masses, positions)
+
+        new_accelerations = []
+        for (acc, acc1) in zip(accelerations, accelerations1):
+            new_acc = []
+            for (a, a1) in zip(acc, acc1):
+                new_acc.append(a + a1)
+            new_accelerations.append(new_acc)
+
+        velocities1 = []
+        for (acc, vel) in zip(new_accelerations, velocities):
+            vel1 = []
+            for (a, v) in zip(acc, vel):
+                vel1.append(v + a*time_step)
+            velocities1.append(vel1)
+        velocities = velocities1
+
+        time += time_step
+
+        if not step % 100:
+            energy, _, _ = compute_energies(masses, positions, velocities)
+            energy_previous = energy
+
+    return energy, energy0
+
+def compute_energies(masses, positions, velocities):
+
+    vel = []
+    for ind, val in enumerate(velocities):
+        vel.append([j**2 for j in val])
+    val = []
+    for i in vel:
+        k = 0
+        for j in i:
+            k += j
+        val.append(k)
+    ke_list = [0.5*i*masses[ind] for ind, i in enumerate(val)]
+    ke = 0.0
+    for i in ke_list:
+        ke += i
+
+    nb_particles = len(masses)
+    pe = 0.0
+    for index_p0 in range(nb_particles - 1):
+        mass0 = masses[index_p0]
+
+        for index_p1 in range(index_p0 + 1, nb_particles):
+
+            mass1 = masses[index_p1]
+            vector = [p0 - p1 for (p0, p1) in zip(positions[index_p0], positions[index_p1])]
+
+            dist = 0
+            for i in vector:
+                dist += i**2
+            distance = math.sqrt(dist)
+
+            pe -= (mass0*mass1) / distance**2
+
+    return ke+pe, ke, pe
+
+if __name__ == "__main__":
+
+    try:
+        time_end = float(sys.argv[2])
+    except IndexError:
+        time_end = 10.0
+
+    time_step = 0.001
+    nb_steps = int(time_end/time_step) + 1
+
+    path_input = sys.argv[1]
+    masses, positions, velocities = load_input_data(path_input)
+
+    masses = masses.tolist()
+    positions = positions.tolist()
+    velocities = velocities.tolist()
+
+    print('time taken', timeit.timeit('loop(time_step, nb_steps, masses, positions, velocities)', globals=globals(), number=50))