diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..7e5a87f
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,14 @@
+OPTIMIZATION ?= 3
+CXXFLAGS     += -O$(OPTIMIZATION)
+LIB          += -lgsl
+
+EXECUTABLE   ?= main
+
+default:
+	$(CXX) $(CPPFLAGS) $(CXXFLAGS) $(INC) loadtxt.cpp main.cpp -o $(EXECUTABLE) $(LIB)
+
+lib:
+	$(CXX) $(CPPFLAGS) $(CXXFLAGS) -fPIC $(INC) loadtxt.cpp main.cpp -shared -o lib$(EXECUTABLE).so $(LIB)
+
+clean:
+	rm -f *.o *.so $(EXECUTABLE)
diff --git a/loadtxt.cpp b/loadtxt.cpp
index 939e0f9..4b60ce5 100644
--- a/loadtxt.cpp
+++ b/loadtxt.cpp
@@ -1,73 +1,69 @@
 #include <algorithm>
 #include <cstdlib>
 #include <fstream>
-#include <vector>
+#include "loadtxt.h"
 
 //TODO if cols is an empty vector, get all columns from the file
 //TODO error checking
 
-class Loadtxt {
-public:
-    Loadtxt(std::string file_name, std::vector<int> cols)
-    {
-        std::sort(cols.begin(), cols.end());
-        n_cols = cols.size();
-        const int tmp_number_of_rows = 16384;
-        int n_rows_alloc = tmp_number_of_rows;
-        buffer = (double*)malloc(n_cols * sizeof(double) * n_rows_alloc);
-        std::ifstream file(file_name);
-        std::string line;
-        int row = -1;
-        while (getline(file, line)) {
-            if (line[line.find_first_not_of(whitespaces)]=='#') continue;
-            if (++row >= n_rows_alloc) {
-                n_rows_alloc += tmp_number_of_rows;
-                buffer = (double*)realloc(buffer, n_cols * sizeof(double) * n_rows_alloc);
-            }
-            line_to_buf(cols, line, buffer + row*n_cols);
+Loadtxt::Loadtxt(std::string file_name, std::vector<int> cols)
+{
+    std::sort(cols.begin(), cols.end());
+    n_cols = cols.size();
+    const int tmp_number_of_rows = 16384;
+    int n_rows_alloc = tmp_number_of_rows;
+    buffer = (double*)malloc(n_cols * sizeof(double) * n_rows_alloc);
+    std::ifstream file(file_name);
+    std::string line;
+    int row = -1;
+    while (getline(file, line)) {
+        if (line[line.find_first_not_of(whitespaces)]=='#') continue;
+        if (++row >= n_rows_alloc) {
+            n_rows_alloc += tmp_number_of_rows;
+            buffer = (double*)realloc(buffer, n_cols * sizeof(double) * n_rows_alloc);
         }
-        file.close();
-        buffer = (double*)realloc(buffer, n_cols * sizeof(double) * (++row));
-        n_rows = row;
+        line_to_buf(cols, line, buffer + row*n_cols);
     }
-    ~Loadtxt()
-    {
-        free(buffer);
-    }
-    std::vector<std::vector<double>> get_cols()
-    {
-        std::vector<std::vector<double>> data(n_cols);
-        for (int col=0; col<n_cols; col++) data[col] = std::vector<double>(n_rows);
-        for (int row=0; row<n_rows; row++) {
-            for (int col=0; col<n_cols; col++) {
-                data[col][row] = buffer[row*n_cols + col];
-            }
+    file.close();
+    buffer = (double*)realloc(buffer, n_cols * sizeof(double) * (++row));
+    n_rows = row;
+}
+
+Loadtxt::~Loadtxt()
+{
+    free(buffer);
+}
+
+std::vector<std::vector<double>> Loadtxt::get_cols()
+{
+    std::vector<std::vector<double>> data(n_cols);
+    for (int col=0; col<n_cols; col++) data[col] = std::vector<double>(n_rows);
+    for (int row=0; row<n_rows; row++) {
+        for (int col=0; col<n_cols; col++) {
+            data[col][row] = buffer[row*n_cols + col];
         }
-        return data;
     }
-private:
-    const char *whitespaces = " \t";
-    void line_to_buf(std::vector<int> cols, std::string line, double *buffer)
-    {
-        int n_cols = cols.size();
-        line = line.substr(line.find_first_not_of(whitespaces));
-        auto pos = line.find_first_of(whitespaces, 1);
-        int col=0, idx=0;
-        std::vector<double> data;
-        while (pos != std::string::npos) {
-            std::string num_as_string;
-            num_as_string = line.substr(0, pos);
-            pos = line.find_first_not_of(whitespaces, pos);
-            line = line.substr(pos, std::string::npos);
-            pos = line.find_first_of(whitespaces, 1);
-            if (col++ == cols[idx]) buffer[idx++] = std::stod(num_as_string);
-        }
-        if (col++ == cols[idx]) buffer[idx++] = std::stod(line);
-        if (idx < n_cols) throw;
+    return data;
+}
+
+void Loadtxt::line_to_buf(std::vector<int> cols, std::string line, double *buffer)
+{
+    int n_cols = cols.size();
+    line = line.substr(line.find_first_not_of(whitespaces));
+    auto pos = line.find_first_of(whitespaces, 1);
+    int col=0, idx=0;
+    std::vector<double> data;
+    while (pos != std::string::npos) {
+        std::string num_as_string;
+        num_as_string = line.substr(0, pos);
+        pos = line.find_first_not_of(whitespaces, pos);
+        line = line.substr(pos, std::string::npos);
+        pos = line.find_first_of(whitespaces, 1);
+        if (col++ == cols[idx]) buffer[idx++] = std::stod(num_as_string);
     }
-    double *buffer;
-    int n_rows, n_cols;
-};
+    if (col++ == cols[idx]) buffer[idx++] = std::stod(line);
+    if (idx < n_cols) throw;
+}
 
 // Below is a deomonstration. The file has multiple columns but we are only
 // interested in the second and fourth. We pass the file name and the column
diff --git a/loadtxt.h b/loadtxt.h
new file mode 100644
index 0000000..02ad1b3
--- /dev/null
+++ b/loadtxt.h
@@ -0,0 +1,14 @@
+#pragma once
+#include <string>
+#include <vector>
+class Loadtxt {
+public:
+    Loadtxt(std::string file_name, std::vector<int> cols);
+    ~Loadtxt();
+    std::vector<std::vector<double>> get_cols();
+private:
+    const char *whitespaces = " \t";
+    void line_to_buf(std::vector<int> cols, std::string line, double *buffer);
+    double *buffer;
+    int n_rows, n_cols;
+};
\ No newline at end of file
diff --git a/main.cpp b/main.cpp
index 2c5ab2b..e2e5624 100644
--- a/main.cpp
+++ b/main.cpp
@@ -1,12 +1,16 @@
-#include <iostream>
-#include <string>
-#include <vector>
+#include <algorithm>
 #include <fstream>
-#include <numeric>
 #include <gsl/gsl_errno.h>
 #include <gsl/gsl_math.h>
 #include <gsl/gsl_odeiv2.h>
 #include <gsl/gsl_spline.h>
+#include <iostream>
+#include <numeric>
+#include <string>
+#include <stdexcept>
+#include <vector>
+
+#include "loadtxt.h"
 
 
 extern "C" const int gsl_success() { return GSL_SUCCESS; } // It's zero, but just for clarity sake.
@@ -22,6 +26,7 @@ public:
         spline = gsl_spline_alloc(gsl_interp_cspline, x.size());
         gsl_spline_init(spline, x.data(), y.data(), x.size()); 
     }
+    Interp() {}
     inline double operator()(double x) const
     {
         return gsl_spline_eval(spline, x, acc);
@@ -54,32 +59,20 @@ class Galaxy {
 public:
     Galaxy(std::string file_name)
     {
-        std::vector<double> t_data, M_halo_data, b_halo_data;
-
-        std::ifstream file(file_name);
-        std::string line;
-        while (std::getline(file, line)) {
-            auto pos = line.find('#');
-            if (pos != std::string::npos) line = line.substr(0, pos);
-            pos = line.find_first_not_of(" \t");
-            if (pos == std::string::npos) continue;
-            double data[3];
-            sscanf(line.c_str(), "%*s %lf %lf %lf", &data[0], &data[1], &data[2]);
-            t_data.push_back(data[0]);
-            M_halo_data.push_back(data[1]);
-            b_halo_data.push_back(data[2]); // note, this is not half-mass radius
-        }
-        interp_M_halo = new Interp(t_data, M_halo_data);
-        interp_b_halo = new Interp(t_data, b_halo_data);
+        auto data = Loadtxt("file.dat", {1, 2, 3}).get_cols();
+        auto& t_data      = data[0];
+        auto& halo_m_data = data[1];
+        auto& halo_b_data = data[2];
+        std::transform(t_data.begin(), t_data.end(), t_data.begin(), [](const double& x){ return x-2.145; });
+        std::transform(halo_b_data.begin(), halo_b_data.end(), halo_b_data.begin(), [](const double& x){ return x*0.7664209365408798; });
+        interp_halo_m = Interp(t_data, halo_m_data);
+        interp_halo_b = Interp(t_data, halo_b_data);
     }
     int func(double t, const double y[], double f[], void *params)
     {
-        double M_halo = (*interp_M_halo)(t);
-        double b_halo = (*interp_b_halo)(t);
-        /*printf("xxxxxxxxx %e, %e msun\n", t, M_halo);
-        printf("xxxxxxxxx %e, %e kpc\n", t, b_halo);
-        exit(0);*/
-        Plummer plummer(M_halo, b_halo);
+        double halo_m = interp_halo_m(t);
+        double halo_b = interp_halo_b(t);
+        Plummer plummer(halo_m, halo_b);
         f[0] = y[3]; // vx -> x'
         f[1] = y[4]; // vy -> y'
         f[2] = y[5]; // vz -> z'
@@ -88,27 +81,25 @@ public:
     }
 
 private:
-    Interp *interp_M_halo;
-    Interp *interp_b_halo;
+    Interp interp_halo_m;
+    Interp interp_halo_b;
 } galaxy("file.dat");
 // Not very nice to have it as a global variable but GSL will have problem otherwise.
 
-int jac(double t, const double y[], double *dfdy, double dfdt[], void *params)
-{
-    return GSL_SUCCESS;
-}
+int jac(double t, const double y[], double *dfdy, double dfdt[], void *params) { return GSL_SUCCESS; }
 
-
-int func(double t, const double y[], double f[], void *params)
+inline int func(double t, const double y[], double f[], void *params)
 {
     return galaxy.func(t, y, f, params);
 }
 
 extern "C"
-int integrate(const double y0[], const double t_max, double y[])
+int integrate(const double y0[], const double t_max, const double step_size, double y[])
 {
-    double t = 2.145;
+    double t = 0;
     constexpr double h = 1./4096.;
+
+    if (step_size/h - (int)(step_size/h) != 0) throw std::runtime_error("step_size must be a multiple of h");
     constexpr double epsabs = 1e-7;
     constexpr double epsrel = 0;
     const gsl_odeiv2_step_type *T = gsl_odeiv2_step_rk8pd;
@@ -118,12 +109,23 @@ int integrate(const double y0[], const double t_max, double y[])
     gsl_odeiv2_system sys = {func, jac, 6, nullptr};
     gsl_odeiv2_driver *d  = gsl_odeiv2_driver_alloc_y_new(&sys, T, h, epsabs, epsrel);
 
+    int step = 0;
+    const int step_max = t_max / step_size;
     std::copy(y0, y0+6, y);
-    int status = gsl_odeiv2_driver_apply(d, &t, t_max, y);
-    return status;
+    for (int step=0; step<step_max; step++) {
+        std::copy(y+step*6, y+(step+1)*6, y+(step+1)*6);
+        int status = gsl_odeiv2_driver_apply(d, &t, (step+1)*step_size, y+(step+1)*6);
+        if (status != GSL_SUCCESS) return status;
+    }
+    return GSL_SUCCESS;
 }
 
 int main()
 {
     std::cout << "bye" << std::endl;
+    double y[12];
+    double y0[] = {80,0,0,0,80,0};
+    for (int i=0; i<30000; i++)
+        integrate(y0, 10, 10, y);
+    return 0;
 }
\ No newline at end of file
diff --git a/plot.py b/plot.py
index e498eb7..dcfa01d 100644
--- a/plot.py
+++ b/plot.py
@@ -4,21 +4,33 @@ import ctypes, subprocess
 recompile   = True
 
 if recompile:
-    p = subprocess.Popen('g++ -shared -o libmain.so -fPIC main.cpp -lgsl'.split())
+    p = subprocess.Popen('g++ -shared -o libmain.so -fPIC loadtxt.cpp main.cpp -lgsl'.split())
     p.wait()
     if p.returncode != 0: raise RuntimeError(p.returncode)
 
 libmain = ctypes.CDLL('./libmain.so')
-def integrate(y0, t_max):
-    y0 = (ctypes.c_double*6)(*y0) 
-    y  = (ctypes.c_double*6)()
-    status = libmain.integrate(y0, ctypes.c_double(t_max), y)
+def integrate(y0, t_max, step_size=None):
+    y0 = (ctypes.c_double*6)(*y0)
+    if step_size is None: step_size = t_max
+    size = int(t_max // step_size) + 1
+    y  = (ctypes.c_double*size*6)()
+    status = libmain.integrate(y0, ctypes.c_double(t_max), ctypes.c_double(step_size), y)
+    y = np.array(y).reshape(size,6)
     return np.array(y), status
 
 gsl_success = libmain.gsl_success()
 
-#t_array = linspace(plot_tmin, plot_tmax, plot_points)
-#r_array = empty_like(t_array)
+from pylab import *
+t_max = 10
+ic = [80,0,0,0,80,0]
+res = integrate(ic, t_max, step_size=1/4096)
+x, y, z, vx, vy, vz = res[0].T
+zzz = x[-1]
+plot(x,y)
 
-res = integrate([10,0,0,0,200,0], 2.245)
-print(res, gsl_success)
\ No newline at end of file
+res = integrate(ic, t_max)
+x, y, z, vx, vy, vz = res[0].T
+plot(x,y,'o')
+print(zzz - x[-1])
+# gca().set_aspect('equal')
+show()
\ No newline at end of file