diff --git a/.gitignore b/.gitignore
index 1c4a1b0..ef08db8 100644
--- a/.gitignore
+++ b/.gitignore
@@ -33,6 +33,11 @@
 
 # Build
 /examples/build/*
+build-*
+cmake-build-*
+
+# IDE
+.idea
 
 # vim temp files
 *.sw*
diff --git a/matplotlibcpp.h b/matplotlibcpp.h
index d95d46a..b12be2c 100644
--- a/matplotlibcpp.h
+++ b/matplotlibcpp.h
@@ -656,6 +656,17 @@ void spy(const std::vector<::std::vector<Numeric>> &x,
 }
 #endif // WITHOUT_NUMPY
 
+class PyAxes {
+public:
+    ~PyAxes() {
+        for (auto &axis: this->axes) {
+            if (axis.second) Py_DECREF(axis.second);
+        }
+    }
+
+    std::map<long, PyObject *> axes;
+};
+
 template <typename Numeric>
 void plot3(const std::vector<Numeric> &x,
                   const std::vector<Numeric> &y,
@@ -723,20 +734,48 @@ void plot3(const std::vector<Numeric> &x,
   }
   if (!fig) throw std::runtime_error("Call to figure() failed.");
 
-  PyObject *gca_kwargs = PyDict_New();
-  PyDict_SetItemString(gca_kwargs, "projection", PyString_FromString("3d"));
-
-  PyObject *gca = PyObject_GetAttrString(fig, "gca");
-  if (!gca) throw std::runtime_error("No gca");
-  Py_INCREF(gca);
-  PyObject *axis = PyObject_Call(
-      gca, detail::_interpreter::get().s_python_empty_tuple, gca_kwargs);
+  static PyAxes axes;
+  PyObject *axis = nullptr;
+  auto it = axes.axes.find(fig_number);
+  if (it != axes.axes.end()) {
+    axis = it->second;
+    Py_INCREF(axis);  // keep it alive for the rest of this function
+  } else {
+    // create new axis
+    PyObject *add_subplot = PyObject_GetAttrString(fig, "add_subplot");
+    if (!add_subplot) throw std::runtime_error("No add_subplot");
+    PyObject *subplot_args = PyTuple_New(1);
+    PyTuple_SetItem(subplot_args, 0, PyLong_FromLong(111));
+    PyObject *subplot_kwargs = PyDict_New();
+    PyDict_SetItemString(subplot_kwargs, "projection", PyString_FromString("3d"));
+
+    axis = PyObject_Call(add_subplot, subplot_args, subplot_kwargs);
+    if (!axis) {
+      PyObject *gca_kwargs = PyDict_New();
+      PyDict_SetItemString(gca_kwargs, "projection", PyString_FromString("3d"));
+
+      PyObject *gca = PyObject_GetAttrString(fig, "gca");
+      if (!gca) throw std::runtime_error("No gca");
+      Py_INCREF(gca);
+      // Fallback for older matplotlib: fig.gca(projection='3d')
+      axis = PyObject_Call(gca, detail::_interpreter::get().s_python_empty_tuple, gca_kwargs);
+      std::cout << "Got old axis: " << axis << std::endl;
+
+      if (!axis) throw std::runtime_error("No axis (both add_subplot and gca failed)");
+
+      Py_INCREF(axis);
+      Py_DECREF(gca);
+      Py_DECREF(gca_kwargs);
+    } else {
+      Py_INCREF(axis);
+    }
 
-  if (!axis) throw std::runtime_error("No axis");
-  Py_INCREF(axis);
+    axes.axes[fig_number] = axis;
 
-  Py_DECREF(gca);
-  Py_DECREF(gca_kwargs);
+    Py_DECREF(add_subplot);
+    Py_DECREF(subplot_args);
+    Py_DECREF(subplot_kwargs);
+  }
 
   PyObject *plot3 = PyObject_GetAttrString(axis, "plot");
   if (!plot3) throw std::runtime_error("No 3D line plot");