diff --git a/examples/2_original_paper.ipynb b/examples/2_original_paper.ipynb
index 4316c7972..53a3ccb31 100644
--- a/examples/2_original_paper.ipynb
+++ b/examples/2_original_paper.ipynb
@@ -44,13 +44,13 @@
     "### Main ideas\n",
     "Suppose we have a set of measurements $x(t)\\in \\mathbb{R}^n$ from some physical system at different points in time $t$. SINDy seeks to represent the time evolution of $x(t)$ in terms of a nonlinear function $f$:\n",
     "\n",
-    "\\begin{equation} \\frac{d}{dt}x(t) = f(x(t)). \\end{equation}\n",
+    "$$ \\frac{d}{dt}x(t) = f(x(t)). $$\n",
     "\n",
     "This equation constitutes a *dynamical system* for the measurements $x(t)$. The vector $x(t)=[x_1(t), x_2(t), \\dots x_n(t)]^\\top$ gives the state of the physical system at time $t$. The function $f(x(t))$ constrains how the system evolves in time.\n",
     "\n",
     "The key idea behind SINDy is that the function $f$ is often **sparse** in the space of an appropriate set of basis functions. For example, the function\n",
     "\n",
-    "\\begin{equation} \\frac{d}{dt}x = f(x) = \\begin{bmatrix} f_1(x)\\\\f_2(x) \\end{bmatrix} = \\begin{bmatrix}1 - x_1 + 3x_1x_2 \\\\ x_2^2 - 5x_1^3 \\end{bmatrix} \\end{equation}\n",
+    "$$ \\frac{d}{dt}x = f(x) = \\begin{bmatrix} f_1(x)\\\\f_2(x) \\end{bmatrix} = \\begin{bmatrix}1 - x_1 + 3x_1x_2 \\\\ x_2^2 - 5x_1^3 \\end{bmatrix} $$\n",
     "\n",
     "is sparse with respect to the set of polynomials of two variables in the sense that if we were to write an expansion of the component functions of $f$ in this basis (e.g. $f_{1}(x) = \\sum_{i=0}^\\infty\\sum_{j=0}^\\infty a_{i,j}x_1^ix_2^j$), only a small number of coefficients ($a_{i,j}$) would be nonzero.\n",
     "\n",
@@ -64,27 +64,27 @@
     "### Approximation problem\n",
     "To apply SINDy in practice one needs a set of measurement data collected at times $t_1, t_2, \\dots, t_n$, and the time derivatives of these measurements (either measured directly or numerically approximated). These data are aggregated into the matrices $X$ and $\\dot X$, respectively\n",
     "\n",
-    "\\begin{equation} X = \\begin{bmatrix} x_1(t_1) & x_2(t_1) & \\dots & x_n(t_1) \\\\ x_1(t_2) & x_2(t_2) & \\dots & x_n(t_2) \\\\ \\vdots & \\vdots & & \\vdots \\\\ x_1(t_m) & x_2(t_m) & \\dots & x_n(t_m) \\\\ \\end{bmatrix}, \\qquad \\dot X = \\begin{bmatrix} \\dot{x_1}(t_1) & \\dot{x_2}(t_1) & \\dots & \\dot{x_n}(t_1) \\\\ \\dot{x_1}(t_2) & \\dot{x_2}(t_2) & \\dots & \\dot{x_n}(t_2) \\\\ \\vdots & \\vdots & & \\vdots \\\\ \\dot{x_1}(t_m) & \\dot{x_2}(t_m) & \\dots & \\dot{x_n}(t_m) \\\\ \\end{bmatrix}. \\end{equation}\n",
+    "$$ X = \\begin{bmatrix} x_1(t_1) & x_2(t_1) & \\dots & x_n(t_1) \\\\ x_1(t_2) & x_2(t_2) & \\dots & x_n(t_2) \\\\ \\vdots & \\vdots & & \\vdots \\\\ x_1(t_m) & x_2(t_m) & \\dots & x_n(t_m) \\\\ \\end{bmatrix}, \\qquad \\dot X = \\begin{bmatrix} \\dot{x_1}(t_1) & \\dot{x_2}(t_1) & \\dots & \\dot{x_n}(t_1) \\\\ \\dot{x_1}(t_2) & \\dot{x_2}(t_2) & \\dots & \\dot{x_n}(t_2) \\\\ \\vdots & \\vdots & & \\vdots \\\\ \\dot{x_1}(t_m) & \\dot{x_2}(t_m) & \\dots & \\dot{x_n}(t_m) \\\\ \\end{bmatrix}. $$\n",
     "\n",
     "Next, one forms a library matrix $\\Theta(X)$ whose columns consist of a chosen set of basis functions applied to the data\n",
     "\n",
-    "\\begin{equation} \\Theta(X) = \\begin{bmatrix} \\mid & \\mid & & \\mid \\\\ \\theta_1(X) & \\theta_2(X) & \\dots & \\theta_\\ell(X) \\\\ \\mid & \\mid & & \\mid  \\end{bmatrix}. \\end{equation}\n",
+    "$$ \\Theta(X) = \\begin{bmatrix} \\mid & \\mid & & \\mid \\\\ \\theta_1(X) & \\theta_2(X) & \\dots & \\theta_\\ell(X) \\\\ \\mid & \\mid & & \\mid  \\end{bmatrix}. $$\n",
     "\n",
     "For example, if $\\theta_1(x), \\theta_2(x), \\dots, \\theta_\\ell(x)$ are monomials ($\\theta_i(x) = x^{i-1}$), then \n",
     "\n",
-    "\\begin{equation} \\theta_3(X) = \\begin{bmatrix} \\mid & \\mid & & \\mid & \\mid & & \\mid \\\\ x_1(t)^2 & x_1(t)x_2(t) & \\dots & x_2(t)^2 & x_2(t)x_3(t) & \\dots & x_n^2(t) \\\\ \\mid & \\mid & & \\mid & \\mid & & \\mid \\end{bmatrix}, \\end{equation}\n",
+    "$$ \\theta_3(X) = \\begin{bmatrix} \\mid & \\mid & & \\mid & \\mid & & \\mid \\\\ x_1(t)^2 & x_1(t)x_2(t) & \\dots & x_2(t)^2 & x_2(t)x_3(t) & \\dots & x_n^2(t) \\\\ \\mid & \\mid & & \\mid & \\mid & & \\mid \\end{bmatrix}, $$\n",
     "\n",
     "where vector products and powers are understood to be element-wise.\n",
     "\n",
     "We seek a set of sparse coefficient vectors (collected into a matrix)\n",
     "\n",
-    "\\begin{equation} \\Xi = \\begin{bmatrix} \\mid & \\mid & & \\mid \\\\ \\xi_1 & \\xi_2 & \\dots & \\xi_n \\\\ \\mid & \\mid & & \\mid  \\end{bmatrix}. \\end{equation}\n",
+    "$$ \\Xi = \\begin{bmatrix} \\mid & \\mid & & \\mid \\\\ \\xi_1 & \\xi_2 & \\dots & \\xi_n \\\\ \\mid & \\mid & & \\mid  \\end{bmatrix}. $$\n",
     "\n",
     "The vector $\\xi_i$ provides the coefficients for a linear combination of basis functions $\\theta_1(x), \\theta_2(x), \\dots, \\theta_\\ell(x)$ representing the $i$th component function of $f$: $f_i(x)$. That is to say, $f_i(x) = \\Theta\\left(x^\\top\\right) \\xi_i$, where $\\Theta\\left(x^\\top\\right)$ is understood to be a row vector consisting of symbolic functions (whereas $\\Theta(X)$ is a matrix whose entries are numerical values).\n",
     "\n",
     "With each of the objects $X$, $\\dot X$, $\\Theta(X)$, and $\\Xi$ being defined, we are ready to write down the approximation problem underlying SINDy:\n",
     "\n",
-    "\\begin{equation} \\dot X \\approx \\Theta(X)\\Xi. \\end{equation}"
+    "$$ \\dot X \\approx \\Theta(X)\\Xi. $$"
    ]
   },
   {
@@ -94,7 +94,7 @@
     "### Structure of PySINDy\n",
     "The submodules of PySINDy are each aligned with one of the terms in the aforementioned approximation equation\n",
     "\n",
-    "\\begin{equation} \\dot X \\approx \\Theta(X)\\Xi. \\end{equation}\n",
+    "$$ \\dot X \\approx \\Theta(X)\\Xi. $$\n",
     "\n",
     " * `pysindy.differentiate` performs numerical differentiation to compute $\\dot X$ from $X$;\n",
     " * `pysindy.feature_library` allows the user to specify a set of library functions and handles the formation of $\\Theta(X)$;\n",
@@ -107,7 +107,7 @@
     "\n",
     "Suppose we have measurements of the position of a particle obeying the following dynamical system at different points in time:\n",
     "\n",
-    "\\begin{equation} \\frac{d}{dt} \\begin{bmatrix} x \\\\ y \\end{bmatrix} = \\begin{bmatrix} -2x \\\\ y \\end{bmatrix} = \\begin{bmatrix} -2 & 0 \\\\ 0 & 1 \\end{bmatrix} \\begin{bmatrix} x \\\\ y \\end{bmatrix} \\end{equation}\n",
+    "$$ \\frac{d}{dt} \\begin{bmatrix} x \\\\ y \\end{bmatrix} = \\begin{bmatrix} -2x \\\\ y \\end{bmatrix} = \\begin{bmatrix} -2 & 0 \\\\ 0 & 1 \\end{bmatrix} \\begin{bmatrix} x \\\\ y \\end{bmatrix} $$\n",
     "\n",
     "Note that this system of differential equations decouples into two differential equations whose solutions are simply $x(t) = x_0e^{-2t}$ and $y(t) = y_0e^t$, where $x_0 = x(0)$ and $y_0=y(0)$ are the initial conditions.\n",
     "\n",
@@ -336,7 +336,7 @@
     "\n",
     "We generate training data by integrating the following linear system of differential equations with initial condtion $(2,0)$.\n",
     "\n",
-    "\\begin{equation} \\frac{d}{dt} \\begin{bmatrix}x \\\\ y\\end{bmatrix} = \\begin{bmatrix} -0.1 & 2 \\\\ -2 & -0.1 \\end{bmatrix} \\begin{bmatrix}x \\\\ y\\end{bmatrix} \\end{equation}"
+    "$$ \\frac{d}{dt} \\begin{bmatrix}x \\\\ y\\end{bmatrix} = \\begin{bmatrix} -0.1 & 2 \\\\ -2 & -0.1 \\end{bmatrix} \\begin{bmatrix}x \\\\ y\\end{bmatrix} $$"
    ]
   },
   {
@@ -464,7 +464,7 @@
     "\n",
     "SINDy can also identify a *nonlinear* two-dimensional damped harmonic oscillator. This example involves cubic dynamics. The governing system used to generate training data is as follows. It is identical to the previous example except the linear terms in the right-hand side have been replaced with cubic ones.\n",
     "\n",
-    "\\begin{equation} \\frac{d}{dt}\\begin{bmatrix}x \\\\ y\\end{bmatrix} = \\begin{bmatrix}-0.1 & 2 \\\\ -2 & -0.1 \\end{bmatrix} \\begin{bmatrix}x^3 \\\\ y^3\\end{bmatrix} \\end{equation}"
+    "$$ \\frac{d}{dt}\\begin{bmatrix}x \\\\ y\\end{bmatrix} = \\begin{bmatrix}-0.1 & 2 \\\\ -2 & -0.1 \\end{bmatrix} \\begin{bmatrix}x^3 \\\\ y^3\\end{bmatrix} $$"
    ]
   },
   {
@@ -576,7 +576,7 @@
     "\n",
     "Our final example of this section shows the application of SINDy to identify of a linear system in three dimensions. The governing equations are given by the following system.\n",
     "\n",
-    "\\begin{equation} \\frac{d}{dt}\\begin{bmatrix}x \\\\ y \\\\ z \\end{bmatrix} = \\begin{bmatrix}-0.1 & -2 & 0 \\\\ 2 & -0.1 & 0 \\\\ 0 & 0 & -0.3 \\end{bmatrix} \\begin{bmatrix}x \\\\ y \\\\ z \\end{bmatrix} \\end{equation}"
+    "$$ \\frac{d}{dt}\\begin{bmatrix}x \\\\ y \\\\ z \\end{bmatrix} = \\begin{bmatrix}-0.1 & -2 & 0 \\\\ 2 & -0.1 & 0 \\\\ 0 & 0 & -0.3 \\end{bmatrix} \\begin{bmatrix}x \\\\ y \\\\ z \\end{bmatrix} $$"
    ]
   },
   {
@@ -710,7 +710,7 @@
     "\n",
     "The Lorenz system serves as an example of a nonlinear ODE whose solutions exhibit chaotic dynamics evolving on a strange attractor. The Lorenz system is given by\n",
     "\n",
-    "\\begin{equation} \\begin{aligned} \\dot x &= \\sigma(y-x)\\\\ \\dot y &= x(\\rho - z) - y \\\\ \\dot z &= xy - \\beta z, \\end{aligned} \\end{equation}\n",
+    "$$ \\begin{aligned} \\dot x &= \\sigma(y-x)\\\\ \\dot y &= x(\\rho - z) - y \\\\ \\dot z &= xy - \\beta z, \\end{aligned} $$\n",
     "\n",
     "with $\\sigma=10$, $\\rho=28$, and $\\beta=\\tfrac{8}{3}$ for this example. We generate our training data starting from the initial condition $(-8, 8, 27)$. Note that for this example we also collect measurements of derivatives of state variables, $(\\dot x, \\dot y, \\dot z)$."
    ]
@@ -1141,7 +1141,7 @@
     "\n",
     "The examples in this section demonstate the ability of SINDy to incorporate parameterized dynamical systems (e.g. normal forms associated with a bifurcation parameter). Note that the parameter, $\\mu$, in these examples is embedded in the dynamics as a variable\n",
     "\n",
-    "\\begin{equation} \\begin{aligned} \\dot{\\mathbf{x}} &= \\mathbf{f}(\\mathbf{x};\\mu) \\\\ \\dot\\mu &= 0. \\end{aligned} \\end{equation}\n",
+    "$$ \\begin{aligned} \\dot{\\mathbf{x}} &= \\mathbf{f}(\\mathbf{x};\\mu) \\\\ \\dot\\mu &= 0. \\end{aligned} $$\n",
     "\n",
     "whose value can the be recovered with SINDy.\n",
     "\n",
@@ -1158,15 +1158,15 @@
     "\n",
     "Discrete time systems require a slightly different treatment than the previous systems. We consider discrete systems of the following form\n",
     "\n",
-    "\\begin{equation} x_{k+1} = f(x_k). \\end{equation}\n",
+    "$$ x_{k+1} = f(x_k). $$\n",
     "\n",
     "To apply SINDy to discrete-time problems, one forms the two data matrices $X_1^{m-1}$ and $X_2^m$, given by\n",
     "\n",
-    "\\begin{equation} X_1^{m-1} = \\begin{bmatrix} - & x_1^\\top & - \\\\ - & x_2^\\top & - \\\\ & \\vdots & \\\\ - & x_{m-1}^\\top & - \\\\ \\end{bmatrix}, \\qquad X_2^{m} = \\begin{bmatrix} - & x_2^\\top & - \\\\ - & x_3^\\top & - \\\\ & \\vdots & \\\\ - & x_{m}^\\top & - \\\\ \\end{bmatrix}. \\end{equation}\n",
+    "$$ X_1^{m-1} = \\begin{bmatrix} - & x_1^\\top & - \\\\ - & x_2^\\top & - \\\\ & \\vdots & \\\\ - & x_{m-1}^\\top & - \\\\ \\end{bmatrix}, \\qquad X_2^{m} = \\begin{bmatrix} - & x_2^\\top & - \\\\ - & x_3^\\top & - \\\\ & \\vdots & \\\\ - & x_{m}^\\top & - \\\\ \\end{bmatrix}. $$\n",
     "\n",
     "The original sparse regression problem is then supplanted by\n",
     "\n",
-    "\\begin{equation} X_2^m \\approx \\Theta\\left(X_1^{m-1}\\right)\\Xi. \\end{equation}"
+    "$$ X_2^m \\approx \\Theta\\left(X_1^{m-1}\\right)\\Xi. $$"
    ]
   },
   {
@@ -1175,7 +1175,7 @@
    "source": [
     "The dynamics of the **logistic map** are described by\n",
     "\n",
-    "\\begin{equation} x_{k+1} = \\mu x_k(1-x_k) + \\eta_k, \\end{equation}\n",
+    "$$ x_{k+1} = \\mu x_k(1-x_k) + \\eta_k, $$\n",
     "\n",
     "where $\\mu$ is the bifurcation parameter and $\\eta_k\\sim N(0, 10^{-3})$ is a (normally distributed) stochastic forcing parameter.\n",
     "\n",
@@ -1396,13 +1396,13 @@
     "\n",
     "We collect noisy data from the Hopf system\n",
     "\n",
-    "\\begin{equation} \\begin{aligned} \\dot x &= \\mu x + \\omega y - Ax(x^2 + y^2) \\\\ \\dot y &= -\\omega x + \\mu y - Ay(x^2 + y^2), \\end{aligned} \\end{equation}\n",
+    "$$ \\begin{aligned} \\dot x &= \\mu x + \\omega y - Ax(x^2 + y^2) \\\\ \\dot y &= -\\omega x + \\mu y - Ay(x^2 + y^2), \\end{aligned} $$\n",
     "\n",
     "with different values of the bifurcation parameter $\\mu$ (fixing $\\omega=A=1$).\n",
     "\n",
     "As in the previous example $\\mu$ is incorporated into the model discovery process by introducing a third equation\n",
     "\n",
-    "\\begin{equation} \\dot \\mu = 0. \\end{equation}\n",
+    "$$ \\dot \\mu = 0. $$\n",
     "\n",
     "If $\\mu<0$ then a (stable) fixed point occurs at $(x,y)=(0,0)$ and if $\\mu>0$ then a (stable) limit cycle occurs on the circle $x^2+y^2 = \\mu$ instead. We test SINDy using two negative values and six positive values of $\\mu$."
    ]