diff --git a/scripts/check_for_memory_leaks.py b/scripts/check_for_memory_leaks.py
new file mode 100755
index 0000000..1f3af61
--- /dev/null
+++ b/scripts/check_for_memory_leaks.py
@@ -0,0 +1,76 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import resource
+import numpy as np
+import ttvfast
+import os
+import matplotlib.pyplot as plt
+
+plt.style.use('ggplot')
+
+def build_args():
+    params = [
+        0.000295994511,
+        0.95573417954,
+        0.00002878248,
+        1.0917340278625494e+01,
+        5.6159310042858110e-02,
+        9.0921164935951211e+01,
+        -1.1729336712101943e-18,
+        1.8094838714599581e+02,
+        -8.7093652691581923e+01,
+        0.00061895914,
+        2.2266898036209028e+01,
+        5.6691301931178648e-02,
+        8.7598285693573246e+01,
+        4.6220554014026838e-01,
+        1.6437004273382669e+00,
+        -1.9584857031843157e+01]
+
+    Time = -1045
+    dt = 0.54
+    Total = 1700
+    n_plan = 2
+    input_flag = 0
+
+    return params, Time, dt, Total, n_plan, input_flag
+
+def run_function(Time, dt, Total, params):
+    # Taken from lweiss test
+    planet1 = ttvfast.models.Planet(*params[2:2+7])
+    planet2 = ttvfast.models.Planet(*params[2+7:])
+
+    gravity, stellar_mass = params[0:2]
+    planets = [planet1, planet2]
+
+    results = ttvfast.ttvfast(planets, stellar_mass, Time, dt, Total)
+
+def rusage_kb():
+    return resource.getrusage(resource.RUSAGE_SELF).ru_maxrss / 1024
+
+def search_for_memory_leaks():
+    args = build_args()
+    Time, dt, Total = args[1:4]
+    dt /= 5.
+    params = args[0]
+
+    iteration, rusages = [], []
+    total = 0
+    for i in range(100):
+        before_usage = rusage_kb()
+        run_function(Time, dt, Total, params)
+        after_usage = rusage_kb()
+        diff_usage = after_usage - before_usage
+        rusages.append(diff_usage + total)
+        total += diff_usage
+        iteration.append(i)
+
+
+    fig, axis = plt.subplots()
+    axis.plot(iteration, rusages)
+    axis.set(title='Total: %d kb' % rusages[-1])
+    plt.show()
+
+if __name__ == '__main__':
+    search_for_memory_leaks()
diff --git a/src/ttvfast_wrapper.c b/src/ttvfast_wrapper.c
index f02771a..81b7458 100644
--- a/src/ttvfast_wrapper.c
+++ b/src/ttvfast_wrapper.c
@@ -11,6 +11,9 @@
 #define PyInteger_FromLong PyInt_FromLong
 #endif
 
+#define STR(x) #x
+#define PRINT_REFCNT(x) (printf("Refcount for %s: %u\n", #x, Py_REFCNT(x)))
+
 void TTVFast(double *params,double dt, double Time, double total,int n_plan,CalcTransit *transit,CalcRV *RV_struct, int nRV, int n_events, int input_flag);
 
 static PyObject *_ttvfast__ttvfast(PyObject *self, PyObject *args);
@@ -94,7 +97,7 @@ init_ttvfast(void)
 
 /* Takes integer and changes the value in place */
 static PyObject *_ttvfast__ttvfast(PyObject *self, PyObject *args) {
-    PyObject *params_obj, *planet_obj, *epoch_obj, *time_obj, *rsky_obj, *vsky_obj, *rv_times_obj, *rv_out_obj;
+    PyObject *params_obj, *planet_obj, *epoch_obj, *time_obj, *rsky_obj, *vsky_obj, *rv_times_obj, *rv_out_obj, *item;
     double dt, Time, total;
     int n_plan, n_events, input_flag, len_rv;
     int i;
@@ -115,11 +118,11 @@ static PyObject *_ttvfast__ttvfast(PyObject *self, PyObject *args) {
     printf("len_rv: %d\n", len_rv);
 #endif
 
+
     /* Get the params list */
     double *params = malloc(sizeof(double) * (2 + n_plan * 7));
-    PyObject *item;
     for (i=0; i<(2 + n_plan * 7); i++) {
-        item = PySequence_GetItem(params_obj, i);
+        item = PyList_GetItem(params_obj, i);
         params[i] = PyFloat_AsDouble(item);
 #ifdef DEBUG
         printf("Params %d: %lf\n", i, params[i]);
@@ -137,8 +140,8 @@ static PyObject *_ttvfast__ttvfast(PyObject *self, PyObject *args) {
     if (len_rv) {
         RV_model = (CalcRV*)calloc(len_rv, sizeof(CalcRV));
         for (int i=0; i<len_rv; i++) {
-            (RV_model+i)->time = PyFloat_AsDouble(
-                    PyList_GetItem(rv_times_obj, i));
+            item = PyList_GetItem(rv_times_obj, i);
+            (RV_model+i)->time = PyFloat_AsDouble(item);
         }
     } else {
         RV_model = NULL;
@@ -157,7 +160,6 @@ static PyObject *_ttvfast__ttvfast(PyObject *self, PyObject *args) {
     time_obj = PyList_New(n_events);
     rsky_obj = PyList_New(n_events);
     vsky_obj = PyList_New(n_events);
-
     rv_out_obj = PyList_New(len_rv);
 
 #ifdef DEBUG
@@ -190,9 +192,21 @@ static PyObject *_ttvfast__ttvfast(PyObject *self, PyObject *args) {
     free(RV_model);
 
     PyObject *positions_out = Py_BuildValue("OOOOO", planet_obj, epoch_obj, time_obj, rsky_obj, vsky_obj);
+
+    Py_XDECREF(planet_obj);
+    Py_XDECREF(epoch_obj);
+    Py_XDECREF(time_obj);
+    Py_XDECREF(rsky_obj);
+    Py_XDECREF(vsky_obj);
+
+    PyObject *retval;
     if (len_rv) {
-        return Py_BuildValue("OO", positions_out, rv_out_obj);
+        retval = Py_BuildValue("OO", positions_out, rv_out_obj);
+        Py_XDECREF(rv_out_obj);
     } else {
-        return Py_BuildValue("OO", positions_out, Py_None);
+        Py_INCREF(Py_None);
+        retval = Py_BuildValue("OO", positions_out, Py_None);
     }
+    Py_XDECREF(positions_out);
+    return retval;
 }