diff --git a/matplotlibcpp.h b/matplotlibcpp.h
index f9f7026b..62b70a1a 100644
--- a/matplotlibcpp.h
+++ b/matplotlibcpp.h
@@ -58,6 +58,7 @@ struct _interpreter {
     PyObject *s_python_function_tight_layout;
     PyObject *s_python_empty_tuple;
     PyObject *s_python_function_stem;
+    PyObject *s_python_function_step;
     PyObject *s_python_function_xkcd;
 
     /* For now, _interpreter is implemented as a singleton since its currently not possible to have
@@ -158,6 +159,7 @@ struct _interpreter {
         s_python_function_errorbar = PyObject_GetAttrString(pymod, "errorbar");
         s_python_function_tight_layout = PyObject_GetAttrString(pymod, "tight_layout");
         s_python_function_stem = PyObject_GetAttrString(pymod, "stem");
+        s_python_function_step = PyObject_GetAttrString(pymod, "step");
         s_python_function_xkcd = PyObject_GetAttrString(pymod, "xkcd");
 
         if(    !s_python_function_show
@@ -187,6 +189,7 @@ struct _interpreter {
             || !s_python_function_errorbar
             || !s_python_function_tight_layout
             || !s_python_function_stem
+            || !s_python_function_step
             || !s_python_function_xkcd
         ) { throw std::runtime_error("Couldn't find required function!"); }
 
@@ -216,6 +219,7 @@ struct _interpreter {
             || !PyFunction_Check(s_python_function_tight_layout)
             || !PyFunction_Check(s_python_function_errorbar)
             || !PyFunction_Check(s_python_function_stem)
+            || !PyFunction_Check(s_python_function_step)
             || !PyFunction_Check(s_python_function_xkcd)
         ) { throw std::runtime_error("Python object is unexpectedly not a PyFunction."); }
 
@@ -337,6 +341,36 @@ bool plot(const std::vector<Numeric> &x, const std::vector<Numeric> &y, const st
     return res;
 }
 
+template<typename Numeric>
+bool step(const std::vector<Numeric> &x, const std::vector<Numeric> &y, const std::map<std::string, std::string>& keywords)
+{
+    assert(x.size() == y.size());
+
+    // using numpy arrays
+    PyObject* xarray = get_array(x);
+    PyObject* yarray = get_array(y);
+
+    // construct positional args
+    PyObject* args = PyTuple_New(2);
+    PyTuple_SetItem(args, 0, xarray);
+    PyTuple_SetItem(args, 1, yarray);
+
+    // construct keyword args
+    PyObject* kwargs = PyDict_New();
+    for(std::map<std::string, std::string>::const_iterator it = keywords.begin(); it != keywords.end(); ++it)
+    {
+        PyDict_SetItemString(kwargs, it->first.c_str(), PyString_FromString(it->second.c_str()));
+    }
+
+    PyObject* res = PyObject_Call(detail::_interpreter::get().s_python_function_step, args, kwargs);
+
+    Py_DECREF(args);
+    Py_DECREF(kwargs);
+    if(res) Py_DECREF(res);
+
+    return res;
+}
+
 template<typename Numeric>
 bool stem(const std::vector<Numeric> &x, const std::vector<Numeric> &y, const std::map<std::string, std::string>& keywords)
 {
@@ -477,6 +511,29 @@ bool plot(const std::vector<NumericX>& x, const std::vector<NumericY>& y, const
     return res;
 }
 
+template<typename NumericX, typename NumericY>
+bool step(const std::vector<NumericX>& x, const std::vector<NumericY>& y, const std::string& s = "")
+{
+    assert(x.size() == y.size());
+
+    PyObject* xarray = get_array(x);
+    PyObject* yarray = get_array(y);
+
+    PyObject* pystring = PyString_FromString(s.c_str());
+
+    PyObject* plot_args = PyTuple_New(3);
+    PyTuple_SetItem(plot_args, 0, xarray);
+    PyTuple_SetItem(plot_args, 1, yarray);
+    PyTuple_SetItem(plot_args, 2, pystring);
+
+    PyObject* res = PyObject_CallObject(detail::_interpreter::get().s_python_function_step, plot_args);
+
+    Py_DECREF(plot_args);
+    if(res) Py_DECREF(res);
+
+    return res;
+}
+
 template<typename NumericX, typename NumericY>
 bool stem(const std::vector<NumericX>& x, const std::vector<NumericY>& y, const std::string& s = "")
 {
@@ -735,6 +792,14 @@ bool plot(const std::vector<Numeric>& y, const std::string& format = "")
     return plot(x,y,format);
 }
 
+template<typename Numeric>
+bool step(const std::vector<Numeric>& y, const std::string& format = "")
+{
+    std::vector<Numeric> x(y.size());
+    for(size_t i=0; i<x.size(); ++i) x.at(i) = i;
+    return step(x,y,format);
+}
+
 template<typename Numeric>
 bool stem(const std::vector<Numeric>& y, const std::string& format = "")
 {