diff --git a/docs/src/examples/beeler_reuter.md b/docs/src/examples/beeler_reuter.md
index 96fc28790..b5ecf7e87 100644
--- a/docs/src/examples/beeler_reuter.md
+++ b/docs/src/examples/beeler_reuter.md
@@ -8,19 +8,23 @@ Note that this tutorial does not use the [higher order IMEX methods built into D
 
 There are hundreds of ionic models that describe cardiac electrical activity in various degrees of detail. Most are based on the classic [Hodgkin-Huxley model](https://en.wikipedia.org/wiki/Hodgkin%E2%80%93Huxley_model) and define the time-evolution of different state variables in the form of nonlinear first-order ODEs. The state vector for these models includes the transmembrane potential, gating variables, and ionic concentrations. The coupling between cells is through the transmembrane potential only and is described as a reaction-diffusion equation, which is a parabolic PDE,
 
-$$\partial V / \partial t = \nabla (D  \nabla V) - \frac {I_\text{ion}} {C_m},$$
+```math
+\frac{\partial V}{\partial t} = \nabla (D  \nabla V) - \frac {I_\text{ion}} {C_m},
+```
 
-where $V$ is the transmembrane potential, $D$ is a diffusion tensor, $I_\text{ion}$ is the sum of the transmembrane currents and is calculated from the ODEs, and $C_m$ is the membrane capacitance, usually assumed to be constant. Here, we model a uniform and isotropic medium. Therefore, the model can be simplified to,
+where ``V`` is the transmembrane potential, ``D`` is a diffusion tensor, ``I_\text{ion}`` is the sum of the transmembrane currents and is calculated from the ODEs, and ``C_m`` is the membrane capacitance, usually assumed to be constant. Here, we model a uniform and isotropic medium. Therefore, the model can be simplified to,
 
-$$\partial V / \partial t = D \Delta{V} - \frac {I_\text{ion}} {C_m},$$
+```math
+\frac{\partial V}{\partial t} = D \nabla^2 V - \frac {I_\text{ion}} {C_m},
+```
 
-where $D$ is now a scalar. By nature, these models have to deal with different time scales and are therefore classified as *stiff*. Commonly, they are solved using the explicit Euler method, typically with a closed form for the integration of the gating variables (the Rush-Larsen method, see below). We can also solve these problems using implicit or semi-implicit PDE solvers (e.g., the [Crank-Nicholson method](https://en.wikipedia.org/wiki/Crank%E2%80%93Nicolson_method) combined with an iterative solver). Higher order explicit methods such as Runge-Kutta and linear multistep methods cannot overcome the stiffness and are not particularly helpful.
+where ``D`` is now a scalar. By nature, these models have to deal with different time scales and are therefore classified as *stiff*. Commonly, they are solved using the explicit Euler method, typically with a closed form for the integration of the gating variables (the Rush-Larsen method, see below). We can also solve these problems using implicit or semi-implicit PDE solvers (e.g., the [Crank-Nicholson method](https://en.wikipedia.org/wiki/Crank%E2%80%93Nicolson_method) combined with an iterative solver). Higher order explicit methods such as Runge-Kutta and linear multistep methods cannot overcome the stiffness and are not particularly helpful.
 
 In this tutorial, we first develop a CPU-only IMEX solver and then show how to move the explicit part to a GPU.
 
 ### The Beeler-Reuter Model
 
-We have chosen the [Beeler-Reuter ventricular ionic model](https://pmc.ncbi.nlm.nih.gov/articles/PMC1283659/) as our example. It is a classic model first described in 1977 and is used as a base for many other ionic models. It has eight state variables, which makes it complicated enough to be interesting without obscuring the main points of the exercise. The eight state variables are: the transmembrane potential ($V$), sodium-channel activation and inactivation gates ($m$ and $h$, similar to the Hodgkin-Huxley model), with an additional slow inactivation gate ($j$), calcium-channel activation and deactivations gates ($d$ and $f$), a time-dependent inward-rectifying potassium current gate ($x_1$), and intracellular calcium concentration ($c$). There are four currents: a sodium current ($i_{Na}$), a calcium current ($i_{Ca}$), and two potassium currents, one time-dependent ($i_{x_1}$) and one background time-independent ($i_{K_1}$).
+We have chosen the [Beeler-Reuter ventricular ionic model](https://pmc.ncbi.nlm.nih.gov/articles/PMC1283659/) as our example. It is a classic model first described in 1977 and is used as a base for many other ionic models. It has eight state variables, which makes it complicated enough to be interesting without obscuring the main points of the exercise. The eight state variables are: the transmembrane potential (``V``), sodium-channel activation and inactivation gates (``m`` and ``h``, similar to the Hodgkin-Huxley model), with an additional slow inactivation gate (``j``), calcium-channel activation and deactivations gates (``d`` and ``f``), a time-dependent inward-rectifying potassium current gate (``x_1``), and intracellular calcium concentration (``c``). There are four currents: a sodium current (``i_\text{Na}``), a calcium current (``i_\text{Ca}``), and two potassium currents, one time-dependent (``i_{x_1}``) and one background time-independent (``i_{K_1}``).
 
 ## CPU-Only Beeler-Reuter Solver
 
@@ -137,27 +141,39 @@ We use an explicit solver for all the state variables except for the transmembra
 
 The [Rush-Larsen](https://ieeexplore.ieee.org/document/4122859/) method replaces the explicit Euler integration for the gating variables with direct integration. The starting point is the general ODE for the gating variables in Hodgkin-Huxley style ODEs,
 
-$$\frac{dg}{dt} = \alpha(V) (1 - g) - \beta(V) g$$
+```math
+\frac{dg}{dt} = (1 - g) \alpha(V) - g \beta(V)
+```
 
-where $g$ is a generic gating variable, ranging from 0 to 1, and $\alpha$ and $\beta$ are reaction rates. This equation can be written as,
+where ``g`` is a generic gating variable, ranging from 0 to 1, and ``α`` and ``β`` are reaction rates. This equation can be written as,
 
-$$\frac{dg}{dt} = (g_{\infty} - g) / \tau_g,$$
+```math
+\frac{dg}{dt} = \frac{g_{\infty} - g}{\tau_g},
+```
 
-where $g_\infty$ and $\tau_g$ are
+where ``g_∞`` and ``τ_g`` are
 
-$$g_{\infty} = \frac{\alpha}{(\alpha + \beta)},$$
+```math
+g_{\infty} = \frac{\alpha}{\alpha + \beta},
+```
 
 and,
 
-$$\tau_g = \frac{1}{(\alpha + \beta)}.$$
+```math
+\tau_g = \frac{1}{\alpha + \beta}.
+```
 
-Assuming that $g_\infty$ and $\tau_g$ are constant for the duration of a single time step ($\Delta{t}$), which is a reasonable assumption for most cardiac models, we can integrate directly to have,
+Assuming that ``g_∞`` and ``τ_g`` are constant for the duration of a single time step (``Δt``), which is a reasonable assumption for most cardiac models, we can integrate directly to have,
 
-$$g(t + \Delta{t}) = g_{\infty} - \left(g_{\infty} - g(\Delta{t})\right)\,e^{-\Delta{t}/\tau_g}.$$
+```math
+g(t + Δt) = g_∞ - \left(g_∞ - g(Δt)\right)\,e^{-Δt/τ_g}.
+```
 
-This is the Rush-Larsen technique. Note that as $\Delta{t} \rightarrow 0$, this equation morphs into the explicit Euler formula,
+This is the Rush-Larsen technique. Note that as ``Δt → 0``, this equation morphs into the explicit Euler formula,
 
-$$g(t + \Delta{t}) = g(t) + \Delta{t}\frac{dg}{dt}.$$
+```math
+g(t + Δt) = g(t) + Δt \frac{dg}{dt}.
+```
 
 `rush_larsen` is a helper function that use the Rush-Larsen method to integrate the gating variables.
 
@@ -169,7 +185,7 @@ $$g(t + \Delta{t}) = g(t) + \Delta{t}\frac{dg}{dt}.$$
 end
 ```
 
-The gating variables are updated as below. The details of how to calculate $\alpha$ and $\beta$ are based on the Beeler-Reuter model and not of direct interest to this tutorial.
+The gating variables are updated as below. The details of how to calculate ``α`` and ``β`` are based on the Beeler-Reuter model and not of direct interest to this tutorial.
 
 ```julia
 function update_M_cpu(g, v, Δt)
@@ -225,8 +241,9 @@ end
 
 ### Implicit Solver
 
-Now, it is time to define the derivative function as an associated function of **BeelerReuterCpu**. We plan to use the CVODE_BDF solver as our implicit portion. Similar to other iterative methods, it calls the deriv function with the same $t$ multiple times. For example, these are consecutive $t$s from a representative run:
+Now, it is time to define the derivative function as an associated function of **BeelerReuterCpu**. We plan to use the `CVODE_BDF` solver as our implicit portion. Similar to other iterative methods, it calls the deriv function with the same ``t`` multiple times. For example, these are consecutive ``t``s from a representative run:
 
+```
 0.86830
 0.86830
 0.85485
@@ -240,8 +257,9 @@ Now, it is time to define the derivative function as an associated function of *
 0.87233
 0.88598
 ...
+```
 
-Here, every time step is called three times. We distinguish between two types of calls to the deriv function. When $t$ changes, the gating variables are updated by calling `update_gates_cpu`:
+Here, every time step is called three times. We distinguish between two types of calls to the deriv function. When ``t`` changes, the gating variables are updated by calling `update_gates_cpu`:
 
 ```julia
 function update_gates_cpu(u, XI, M, H, J, D, F, C, Δt)
@@ -265,7 +283,7 @@ function update_gates_cpu(u, XI, M, H, J, D, F, C, Δt)
 end
 ```
 
-On the other hand, du is updated at each time step, since it is independent of $\Delta{t}$.
+On the other hand, du is updated at each time step, since it is independent of ``Δt``.
 
 ```julia
 # iK1 is the inward-rectifying potassium current
@@ -367,7 +385,7 @@ deriv_cpu = BeelerReuterCpu(u0, 1.0);
 prob = DE.ODEProblem(deriv_cpu, u0, (0.0, 50.0));
 ```
 
-For stiff reaction-diffusion equations, CVODE_BDF from Sundial library is an excellent solver.
+For stiff reaction-diffusion equations, `CVODE_BDF` from Sundial library is an excellent solver.
 
 ```julia
 @time sol = DE.solve(prob, Sundials.CVODE_BDF(linear_solver = :GMRES), saveat = 100.0);
@@ -405,7 +423,7 @@ Some libraries, such as [ArrayFire](https://github.com/arrayfire/arrayfire), hid
 
 The key to fast CUDA programs is to minimize CPU/GPU memory transfers and global memory accesses. The implicit solver is currently CPU only, but it only requires access to the transmembrane potential. The rest of state variables reside on the GPU memory.
 
-We modify ``BeelerReuterCpu`` into ``BeelerReuterGpu`` by defining the state variables as *CuArray*s instead of standard Julia *Array*s. The name of each variable defined on the GPU is prefixed by *d_* for clarity. Note that $\Delta{v}$ is a temporary storage for the Laplacian and stays on the CPU side.
+We modify `BeelerReuterCpu` into `BeelerReuterGpu` by defining the state variables as *CuArray*s instead of standard Julia *Array*s. The name of each variable defined on the GPU is prefixed by *d_* for clarity. Note that ``Δv`` is a temporary storage for the Laplacian and stays on the CPU side.
 
 ```julia
 import CUDA
@@ -556,7 +574,7 @@ A CUDA program does not directly deal with GPCs and SMs. The logical view of a C
 
 Each thread can find its logical coordinate by using few pre-defined indexing variables (*threadIdx*, *blockIdx*, *blockDim* and *gridDim*) in C/C++ and the corresponding functions (e.g., `threadIdx()`) in Julia. Their variables and functions are defined automatically for each thread and may return a different value depending on the calling thread. The return value of these functions is a 1-, 2-, or 3-dimensional structure whose elements can be accessed as `.x`, `.y`, and `.z` (for a 1-dimensional case, `.x` reports the actual index and `.y` and `.z` simply return 1). For example, if we deploy a kernel in 128 blocks and with 256 threads per block, each thread will see
 
-```
+```c
     gridDim.x = 128;
     blockDim=256;
 ```
@@ -565,13 +583,13 @@ while `blockIdx.x` ranges from 0 to 127 in C/C++ and 1 to 128 in Julia. Similarl
 
 A C/C++ thread can calculate its index as
 
-```
+```c
     int idx = blockDim.x * blockIdx.x + threadIdx.x;
 ```
 
 In Julia, we have to take into account base 1. Therefore, we use the following formula
 
-```
+```julia
     idx = (blockIdx().x-UInt32(1)) * blockDim().x + threadIdx().x
 ```
 
diff --git a/docs/src/examples/classical_physics.md b/docs/src/examples/classical_physics.md
index dfe00ee01..4170ceb1d 100644
--- a/docs/src/examples/classical_physics.md
+++ b/docs/src/examples/classical_physics.md
@@ -6,9 +6,11 @@ If you're getting some cold feet to jump in to DiffEq land, here are some handcr
 
 #### Radioactive Decay of Carbon-14
 
-$$f(t,u) = \frac{du}{dt}$$
+```math
+\frac{du}{dt} = f(t,u) = - λ u
+```
 
-The Radioactive decay problem is the first order linear ODE problem of an exponential with a negative coefficient, which represents the half-life of the process in question. Should the coefficient be positive, this would represent a population growth equation.
+The Radioactive decay problem is the first order linear ODE problem of an exponential with a negative coefficient, which represents the half-life of the process in question. Should the coefficient be positive, this would represent a population growth equation. ``λ`` is known as the decay rate, and can be related to the half-life as ``λ = \ln(2)/t_{1/2}``.
 
 ```@example physics
 import OrdinaryDiffEq as ODE, Plots
@@ -41,21 +43,27 @@ Plots.plot!(sol.t, t -> 2^(-t / t½), lw = 3, ls = :dash, label = "Analytical So
 
 Another classical example is the harmonic oscillator, given by:
 
-$$\ddot{x} + \omega^2 x = 0$$
+```math
+\ddot{x} + \omega^2 x = 0
+```
 
 with the known analytical solution
 
-$$\begin{align*}
+```math
+\begin{align*}
 x(t) &= A\cos(\omega t - \phi) \\
 v(t) &= -A\omega\sin(\omega t - \phi),
-\end{align*}$$
+\end{align*}
+```
 
 where
 
-$$A = \sqrt{c_1 + c_2} \qquad\text{and}\qquad \tan \phi = \frac{c_2}{c_1}$$
+```math
+A = \sqrt{c_1 + c_2} \qquad\text{and}\qquad \tan \phi = \frac{c_2}{c_1}
+```
 
-with $c_1, c_2$ constants determined by the initial conditions such that
-$c_1$ is the initial position and $\omega c_2$ is the initial velocity.
+with ``c_1``, ``c_2`` constants determined by the initial conditions such that
+``c_1`` is the initial position and ``\omega c_2`` is the initial velocity.
 
 Instead of transforming this to a system of ODEs to solve with `ODEProblem`,
 we can use `SecondOrderODEProblem` as follows.
@@ -98,22 +106,28 @@ Thus, if we want the first series to be `x`, we have to flip the order with `var
 
 #### Simple Pendulum
 
-We will start by solving the pendulum problem. In the physics class, we often solve this problem by small angle approximation, i.e. $ sin(\theta) \approx \theta$, because otherwise, we get an elliptic integral which doesn't have an analytic solution. The linearized form is
+We will start by solving the pendulum problem. In the physics class, we often solve this problem by small angle approximation, i.e. ``\sin(θ) \approx θ``, because otherwise, we get an elliptic integral which doesn't have an analytic solution. The linearized form is
 
-$$\ddot{\theta} + \frac{g}{L}{\theta} = 0$$
+```math
+\ddot{θ} + \frac{g}{L} θ = 0
+```
 
 But we have numerical ODE solvers! Why not solve the *real* pendulum?
 
-$$\ddot{\theta} + \frac{g}{L}{\sin(\theta)} = 0$$
+```math
+\ddot{θ} + \frac{g}{L} \sin(θ) = 0
+```
 
 Notice that now we have a second order ODE.
 In order to use the same method as above, we need to transform it into a system
-of first order ODEs by employing the notation $d\theta = \dot{\theta}$.
+of first order ODEs by employing the notation ``ω(t) = \dot{θ}``.
 
-$$\begin{align*}
-&\dot{\theta} = d{\theta} \\
-&\dot{d\theta} = - \frac{g}{L}{\sin(\theta)}
-\end{align*}$$
+```math
+\begin{align*}
+\dot{θ} &= ω \\
+\dot{ω} &= - \frac{g}{L} \sin(θ)
+\end{align*}
+```
 
 ```@example physics
 # Simple Pendulum Problem
@@ -129,9 +143,8 @@ tspan = (0.0, 6.3)
 
 #Define the problem
 function simplependulum(du, u, p, t)
-    θ = u[1]
-    dθ = u[2]
-    du[1] = dθ
+    θ, ω = u
+    du[1] = ω
     du[2] = -(g / L) * sin(θ)
 end
 
@@ -167,16 +180,17 @@ Plots.plot(p, xlims = (-9, 9))
 A more complicated example is given by the double pendulum. The equations governing
 its motion are given by the following (taken from this [Stack Overflow question](https://mathematica.stackexchange.com/questions/40122/help-to-plot-poincar%C3%A9-section-for-double-pendulum))
 
-$$\frac{d}{dt}
-\begin{pmatrix}
-\alpha \\ l_\alpha \\ \beta \\ l_\beta
-\end{pmatrix}=
+```math
+\frac{d}{dt}
+\begin{pmatrix} \alpha \\ l_\alpha \\ \beta \\ l_\beta \end{pmatrix}
+=
 \begin{pmatrix}
 2\frac{l_\alpha - (1+\cos\beta)l_\beta}{3-\cos 2\beta} \\
 -2\sin\alpha - \sin(\alpha + \beta) \\
 2\frac{-(1+\cos\beta)l_\alpha + (3+2\cos\beta)l_\beta}{3-\cos2\beta}\\
 -\sin(\alpha+\beta) - 2\sin(\beta)\frac{(l_\alpha-l_\beta)l_\beta}{3-\cos2\beta} + 2\sin(2\beta)\frac{l_\alpha^2-2(1+\cos\beta)l_\alpha l_\beta + (3+2\cos\beta)l_\beta^2}{(3-\cos2\beta)^2}
-\end{pmatrix}$$
+\end{pmatrix}
+```
 
 ```@example physics
 #Double Pendulum Problem
@@ -197,7 +211,7 @@ function polar2cart(sol; dt = 0.02, l1 = L₁, l2 = L₂, vars = (2, 4))
     x1 = l1 * sin.(p1)
     y1 = l1 * -cos.(p1)
     (u, (x1 + l2 * sin.(p2),
-        y1 - l2 * cos.(p2)))
+         y1 - l2 * cos.(p2)))
 end
 
 #Define the Problem
@@ -231,7 +245,7 @@ Instead of looking at the full phase space, we can look at Poincaré sections,
 which are sections through a higher-dimensional phase space diagram.
 This helps to understand the dynamics of interactions and is wonderfully pretty.
 
-The Poincaré section in this is given by the collection of $(β,l_β)$ when $α=0$ and $\frac{dα}{dt}>0$.
+The Poincaré section in this is given by the collection of ``(β,l_β)`` when ``α=0`` and ``\frac{dα}{dt}>0``.
 
 ```@example physics
 #Constants and setup
@@ -241,10 +255,7 @@ tspan2 = (0.0, 500.0)
 
 #Define the problem
 function double_pendulum_hamiltonian(udot, u, p, t)
-    α = u[1]
-    lα = u[2]
-    β = u[3]
-    lβ = u[4]
+    α, lα, β, lβ = u
     udot .= [2(lα - (1 + cos(β))lβ) / (3 - cos(2β)),
         -2sin(α) - sin(α + β),
         2(-(1 + cos(β))lα + (3 + 2cos(β))lβ) / (3 - cos(2β)),
@@ -284,22 +295,30 @@ Plots.plot(p, xlabel = "\\beta", ylabel = "l_\\beta", ylims = (0, 0.03))
 
 The Hénon-Heiles potential occurs when non-linear motion of a star around a galactic center, with the motion restricted to a plane.
 
-$$\begin{align}
+```math
+\begin{align*}
 \frac{d^2x}{dt^2}&=-\frac{\partial V}{\partial x}\\
 \frac{d^2y}{dt^2}&=-\frac{\partial V}{\partial y}
-\end{align}$$
+\end{align*}
+```
 
 where
 
-$$V(x,y)={\frac {1}{2}}(x^{2}+y^{2})+\lambda \left(x^{2}y-{\frac {y^{3}}{3}}\right).$$
+```math
+V(x,y) = \frac {1}{2} (x^2+y^2) + λ \left(x^2 y - \frac{y^3}{3}\right).
+```
 
-We pick $\lambda=1$ in this case, so
+We pick ``λ=1`` in this case, so
 
-$$V(x,y) = \frac{1}{2}(x^2+y^2+2x^2y-\frac{2}{3}y^3).$$
+```math
+V(x,y) = \frac{1}{2} \left(x^2 + y^2 + 2 x^2 y - \frac{2}{3} y^3\right).
+```
 
 Then the total energy of the system can be expressed by
 
-$$E = T+V = V(x,y)+\frac{1}{2}(\dot{x}^2+\dot{y}^2).$$
+```math
+E = T+V = V(x,y) + \frac{1}{2} \left(\dot{x}^2+\dot{y}^2\right).
+```
 
 The total energy should conserve as this system evolves.
 
@@ -317,10 +336,7 @@ E(x, y, dx, dy) = V(x, y) + 1 // 2 * (dx^2 + dy^2);
 
 #Define the function
 function Hénon_Heiles(du, u, p, t)
-    x = u[1]
-    y = u[2]
-    dx = u[3]
-    dy = u[4]
+    x, y, dx, dy = u
     du[1] = dx
     du[2] = dy
     du[3] = -x - 2x * y
diff --git a/docs/src/examples/diffusion_implicit_heat_equation.md b/docs/src/examples/diffusion_implicit_heat_equation.md
index 237603124..af6f2acff 100644
--- a/docs/src/examples/diffusion_implicit_heat_equation.md
+++ b/docs/src/examples/diffusion_implicit_heat_equation.md
@@ -3,10 +3,10 @@
 In this tutorial, we'll be solving the heat equation:
 
 ```math
-∂_t T = α ∇²(T) + β \sin(γ z)
+∂_t T = α ∇²T + β \sin(γ z)
 ```
 
-with boundary conditions: ``∇T(z=a) = ∇T_{bottom}, T(z=b) = T_{top}``. We'll solve these equations numerically using Finite Difference Method on cell faces. The same exercise could easily be done on cell centers.
+with boundary conditions: ``∇T(z=a) = ∇T_\text{bottom}`` and ``T(z=b) = T_\text{top}``. We'll solve these equations numerically using Finite Difference Method on cell faces. The same exercise could easily be done on cell centers.
 
 ## Code loading and parameters
 
@@ -44,23 +44,29 @@ zf = range(a, b, length = n + 1);   # coordinates on cell faces
 Here, we'll derive the analytic solution:
 
 ```math
-\frac{∂²T}{∂²z} = -\frac{S(z)}{α} = -\frac{β \sin(γ z)}{α} \\
-\frac{∂T}{∂z} = \frac{β \cos(γ z)}{γ α}+c_1 \\
-T(z) = \frac{β \sin(γ z)}{γ^2 α}+c_1 z+c_2, \qquad \text{(generic solution)}
+\begin{align*}
+\frac{∂²T}{∂²z} &= -\frac{S(z)}{α} = -\frac{β \sin(γ z)}{α} \\
+\frac{∂T}{∂z} &= \frac{β \cos(γ z)}{γ α}+c_1 \\
+T(z) &= \frac{β \sin(γ z)}{γ^2 α}+c_1 z+c_2, \qquad \text{(generic solution)}
+\end{align*}
 ```
 
 Apply bottom boundary condition:
 
 ```math
-\frac{∂T}{∂z}(a) = \frac{β \cos(γ a)}{γ α}+c_1 = ∇T_{bottom} \\
-c_1 = ∇T_{bottom}-\frac{β \cos(γ a)}{γ α}
+\begin{align*}
+\frac{∂T}{∂z}(a) &= \frac{β \cos(γ a)}{γ α} + c_1 = ∇T_{bottom} \\
+c_1 &= ∇T_\text{bottom} - \frac{β \cos(γ a)}{γ α}
+\end{align*}
 ```
 
 Apply top boundary condition:
 
 ```math
-T(b) = \frac{β \sin(γ b)}{γ^2 α}+c_1 b+c_2 = T_{top} \\
-c_2 = T_{top}-\left(\frac{β \sin(γ b)}{γ^2 α}+c_1 b\right)
+\begin{align*}
+T(b) &= \frac{β \sin(γ b)}{γ^2 α} + c_1 b+c_2 = T_\text{top} \\
+c_2 &= T_\text{top} - \left(\frac{β \sin(γ b)}{γ^2 α}+c_1 b\right)
+\end{align*}
 ```
 
 And now let's define this in a Julia function:
@@ -78,10 +84,12 @@ end
 Here, we'll derive the matrix form of the temporal discretization we wish to use (diffusion-implicit and explicit Euler):
 
 ```math
-∂_t T = α ∇²T + S \\
-(T^{n+1}-T^n) = Δt (α  ∇²T^{n+1} + S) \\
-(T^{n+1} - Δt α ∇²T^{n+1}) = T^n + Δt S \\
-(I - Δt α ∇²) T^{n+1} = T^n + Δt S
+\begin{align*}
+∂_t T &= α ∇²T + S \\
+(T^{n+1}-T^n) &= Δt (α  ∇²T^{n+1} + S) \\
+(T^{n+1} - Δt α ∇²T^{n+1}) &= T^n + Δt S \\
+(I - Δt α ∇²) T^{n+1} &= T^n + Δt S
+\end{align*}
 ```
 
 Note that, since the ``∇²`` reaches to boundary points, we'll need to modify the stencils to account for boundary conditions.
@@ -91,8 +99,10 @@ Note that, since the ``∇²`` reaches to boundary points, we'll need to modify
 For the interior domain, a central and second-order finite difference stencil is simply:
 
 ```math
-∇²f = \frac{f_{i-1} -2f_i + f_{i+1}}{Δz²}, \qquad \text{or} \\
-∇² = \left[\frac{1}{Δz²}, \frac{-2}{Δz²}, \frac{1}{Δz²}\right] \\
+\begin{align*}
+∇²f &= \frac{f_{i-1} -2f_i + f_{i+1}}{Δz²}, \qquad \text{or} \\
+∇² &= \left[\frac{1}{Δz²}, \frac{-2}{Δz²}, \frac{1}{Δz²}\right] \\
+\end{align*}
 ```
 
 At the boundaries, we need to modify the stencil to account for Dirichlet and Neumann BCs. Using the following index denotation:
@@ -104,18 +114,22 @@ At the boundaries, we need to modify the stencil to account for Dirichlet and Ne
 the Dirichlet boundary stencil & source:
 
 ```math
-∂_t T = α \frac{T[i-1]+T[b]-2 T[i]}{Δz²} + S \\
-∂_t T = α \frac{T[i-1]-2 T[i]}{Δz²} + S + α \frac{T[b]}{Δz²}
+\begin{align*}
+∂_t T &= α \frac{T[i-1]+T[b]-2 T[i]}{Δz²} + S \\
+∂_t T &= α \frac{T[i-1]-2 T[i]}{Δz²} + S + α \frac{T[b]}{Δz²}
+\end{align*}
 ```
 
 and Neumann boundary stencil & source:
 
 ```math
-∇T_{bottom} n̂ = \frac{T[g] - T[i]}{2Δz}, \qquad    n̂ = [-1,1] ∈ [z_{min},z_{max}] \\
-T[i] + 2 Δz ∇T_{bottom} n̂ = T[g] \\
-∂_t T = α \frac{\frac{(T[i] + 2 Δz ∇T_{bottom} n̂) - T[b]}{Δz} - \frac{T[b] - T[i]}{Δz}}{Δz} + S \\
-∂_t T = α \frac{\frac{T[i] - T[b]}{Δz} - \frac{T[b] - T[i]}{Δz}}{Δz} + S + α 2 Δz \frac{∇T_{bottom}}{Δz²} \\
-∂_t T = α \frac{2 T[i] - 2 T[b]}{Δz²} + S + 2α \frac{∇T_{bottom} n̂}{Δz}
+\begin{align*}
+∇T_\text{bottom} n̂ = \frac{T[g] - T[i]}{2Δz}, \qquad    n̂ = [-1,1] ∈ [z_\text{min},z_\text{max}] \\
+T[i] + 2 Δz ∇T_\text{bottom} n̂ = T[g] \\
+∂_t T &= α \frac{\frac{(T[i] + 2 Δz ∇T_\text{bottom} n̂) - T[b]}{Δz} - \frac{T[b] - T[i]}{Δz}}{Δz} + S \\
+∂_t T &= α \frac{\frac{T[i] - T[b]}{Δz} - \frac{T[b] - T[i]}{Δz}}{Δz} + S + α 2 Δz \frac{∇T_\text{bottom}}{Δz²} \\
+∂_t T &= α \frac{2 T[i] - 2 T[b]}{Δz²} + S + 2α \frac{∇T_\text{bottom} n̂}{Δz}
+\end{align*}
 ```
 
 ## Define the discrete diffusion operator
diff --git a/docs/src/examples/kepler_problem.md b/docs/src/examples/kepler_problem.md
index edefaf7b9..9b1c56265 100644
--- a/docs/src/examples/kepler_problem.md
+++ b/docs/src/examples/kepler_problem.md
@@ -19,8 +19,9 @@ Also, we know that
 ```
 
 ```@example kepler
-import OrdinaryDiffEq as ODE, LinearAlgebra, ForwardDiff, NonlinearSolve as NLS, Plots
-H(q, p) = LinearAlgebra.norm(p)^2 / 2 - inv(LinearAlgebra.norm(q))
+import OrdinaryDiffEq as ODE, ForwardDiff, Plots
+import LinearAlgebra: norm
+H(q, p) = norm(p)^2 / 2 - inv(norm(q))
 L(q, p) = q[1] * p[2] - p[1] * q[2]
 
 pdot(dp, p, q, params, t) = ForwardDiff.gradient!(dp, q -> -H(q, p), q)
@@ -35,10 +36,6 @@ prob = ODE.DynamicalODEProblem(pdot, qdot, initial_velocity, initial_position, t
 sol = ODE.solve(prob, ODE.KahanLi6(), dt = 1 // 10);
 ```
 
-!!! note
-    
-    Note that NonlinearSolve.jl is required to be imported for ManifoldProjection
-
 Let's plot the orbit and check the energy and angular momentum variation. We know that energy and angular momentum should be constant, and they are also called first integrals.
 
 ```@example kepler
@@ -103,8 +100,12 @@ Both Runge-Kutta-Nyström and Runge-Kutta integrator do not conserve energy nor
 
 In this example, we know that energy and angular momentum should be conserved. We can achieve this through manifold projection. As the name implies, it is a procedure to project the ODE solution to a manifold. Let's start with a base case, where manifold projection isn't being used.
 
+!!! note
+    
+    Note that NonlinearSolve.jl is required to be imported for ManifoldProjection
+
 ```@example kepler
-import DiffEqCallbacks as CB
+import DiffEqCallbacks as CB, NonlinearSolve as NLS
 
 function plot_orbit2(sol)
     Plots.plot(sol, vars = (1, 2), lab = "Orbit", title = "Kepler Problem Solution")
diff --git a/docs/src/examples/min_and_max.md b/docs/src/examples/min_and_max.md
index 7ff536bd2..6209e2ce9 100644
--- a/docs/src/examples/min_and_max.md
+++ b/docs/src/examples/min_and_max.md
@@ -14,10 +14,7 @@ tspan = (0.0, 100.0)
 
 #Define the problem
 function double_pendulum_hamiltonian(udot, u, p, t)
-    α = u[1]
-    lα = u[2]
-    β = u[3]
-    lβ = u[4]
+    α, lα, β, lβ = u
     udot .= [2(lα - (1 + cos(β))lβ) / (3 - cos(2β)),
         -2sin(α) - sin(α + β),
         2(-(1 + cos(β))lα + (3 + 2cos(β))lβ) / (3 - cos(2β)),
diff --git a/docs/src/examples/outer_solar_system.md b/docs/src/examples/outer_solar_system.md
index 91e5611ee..3d1170974 100644
--- a/docs/src/examples/outer_solar_system.md
+++ b/docs/src/examples/outer_solar_system.md
@@ -66,7 +66,7 @@ potential = -G *
 `NBodyProblem` constructs a second order ODE problem under the hood. We know that a Hamiltonian system has the form of
 
 ```math
-\dot{p} = -H_{q}(p,q), \quad \dot{q} = H_{p}(p,q)
+\dot{p} = -\frac{\partial H}{\partial q}, \quad \dot{q} = \frac{\partial H}{\partial p}
 ```
 
 For an N-body system, we can simplify this as:
diff --git a/docs/src/examples/spiking_neural_systems.md b/docs/src/examples/spiking_neural_systems.md
index 944b9ae12..ba04bca8b 100644
--- a/docs/src/examples/spiking_neural_systems.md
+++ b/docs/src/examples/spiking_neural_systems.md
@@ -185,7 +185,7 @@ function HH!(du, u, p, t)
     gK, gNa, gL, EK, ENa, EL, C, I = p
     v, n, m, h = u
 
-    du[1] = (-(gK * (n^4.0) * (v - EK)) - (gNa * (m^3.0) * h * (v - ENa)) -
+    du[1] = (-(gK * n^4 * (v - EK)) - (gNa * m^3 * h * (v - ENa)) -
              (gL * (v - EL)) + I) / C
     du[2] = (alpha_n(v) * (1.0 - n)) - (beta_n(v) * n)
     du[3] = (alpha_m(v) * (1.0 - m)) - (beta_m(v) * m)
diff --git a/docs/src/getting_started.md b/docs/src/getting_started.md
index f49138ed2..4150fda54 100644
--- a/docs/src/getting_started.md
+++ b/docs/src/getting_started.md
@@ -375,8 +375,8 @@ Parameterized functions can also be used for building **nonhomogeneous ordinary
 ```math
 \begin{aligned}
 \frac{\mathrm{d}\theta(t)}{\mathrm{d}t} &= \omega(t)\\
-\frac{\mathrm{d}\omega(t)}{\mathrm{d}t} &= - \frac{3}{2}\frac{g}{l}\sin\theta(t) + \frac{3}{ml^2}M(t)
-\end{aligned},
+\frac{\mathrm{d}\omega(t)}{\mathrm{d}t} &= - \frac{3}{2}\frac{g}{l}\sin\theta(t) + \frac{3}{ml^2}M(t),
+\end{aligned}
 ```
 
 where `θ` and `ω` are the angular deviation of the pendulum from the vertical (hanging) orientation and the angular rate, respectively, `M` is an external torque (developed, say, by a wind or a motor), and finally, `g` stands for gravitational acceleration.
@@ -433,7 +433,7 @@ f(u, p, t) = A * u
 prob = DE.ODEProblem(f, u0, tspan)
 ```
 
-Here our ODE is on a 4x2 matrix, and the ODE is the linear system defined by
+Here our ODE is on a 4×2 matrix, and the ODE is the linear system defined by
 multiplication by `A`. To solve the ODE, we do the same steps
 as before.
 
diff --git a/docs/src/tutorials/advanced_ode_example.md b/docs/src/tutorials/advanced_ode_example.md
index f410f9813..b49564858 100644
--- a/docs/src/tutorials/advanced_ode_example.md
+++ b/docs/src/tutorials/advanced_ode_example.md
@@ -25,10 +25,10 @@ differential equation (BRUSS).
 The Brusselator PDE is defined on a unit square periodic domain as follows:
 
 ```math
-\begin{align}
+\begin{align*}
 \frac{\partial U}{\partial t} &= 1 + U^2V - 4.4U + \alpha \nabla^2 U + f(x, y, t),\\
 \frac{\partial V}{\partial t} &= 3.4U - U^2V + \alpha \nabla^2 V,
-\end{align}
+\end{align*}
 ```
 
 where
@@ -36,26 +36,26 @@ where
 ```math
 f(x, y, t) = \begin{cases}
 5 & \quad \text{if } (x-0.3)^2+(y-0.6)^2 ≤ 0.1^2 \text{ and } t ≥ 1.1\\
-0 & \quad \text{else}
-\end{cases}, \mathrm{and}
+0 & \quad \text{else},
+\end{cases}
 ```
 
-$\nabla^2 = \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}$ is the two dimensional Laplacian operator. The above equations are to be solved for a time interval $t \in [0, 11.5]$ subject to the initial conditions
+and ``\nabla^2 = \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}`` is the two dimensional Laplacian operator. The above equations are to be solved for a time interval ``t \in [0, 11.5]`` subject to the initial conditions
 
 ```math
-\begin{align}
+\begin{align*}
 U(x, y, 0) &= 22\cdot (y(1-y))^{3/2} \\
-V(x, y, 0) &= 27\cdot (x(1-x))^{3/2}
-\end{align},
+V(x, y, 0) &= 27\cdot (x(1-x))^{3/2},
+\end{align*}
 ```
 
 and the periodic boundary conditions
 
 ```math
-\begin{align}
+\begin{align*}
 U(x+1,y,t) &= U(x,y,t) \\
 V(x,y+1,t) &= V(x,y,t).
-\end{align}
+\end{align*}
 ```
 
 To solve this PDE, we will discretize it into a system of ODEs with the finite
@@ -295,10 +295,8 @@ function algebraicmultigrid2(W, du, u, p, t, newW, Plprev, Prprev, solverdata)
     if newW === nothing || newW
         A = convert(AbstractMatrix, W)
         Pl = AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.ruge_stuben(A,
-            presmoother = AlgebraicMultigrid.Jacobi(rand(size(A,
-                1))),
-            postsmoother = AlgebraicMultigrid.Jacobi(rand(size(A,
-                1)))))
+            presmoother = AlgebraicMultigrid.Jacobi(rand(size(A, 1))),
+            postsmoother = AlgebraicMultigrid.Jacobi(rand(size(A, 1)))))
     else
         Pl = Plprev
     end
@@ -415,10 +413,8 @@ And similarly for algebraic multigrid:
 
 ```julia
 prectmp2 = AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.ruge_stuben(W,
-    presmoother = AlgebraicMultigrid.Jacobi(rand(size(W,
-        1))),
-    postsmoother = AlgebraicMultigrid.Jacobi(rand(size(W,
-        1)))))
+    presmoother = AlgebraicMultigrid.Jacobi(rand(size(W, 1))),
+    postsmoother = AlgebraicMultigrid.Jacobi(rand(size(W, 1)))))
 const preccache2 = Ref(prectmp2)
 function psetupamg(p, t, u, du, jok, jcurPtr, gamma)
     if jok
@@ -433,10 +429,8 @@ function psetupamg(p, t, u, du, jok, jcurPtr, gamma)
         # Build preconditioner on W
         preccache2[] = AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.ruge_stuben(
             W,
-            presmoother = AlgebraicMultigrid.Jacobi(rand(size(W,
-                1))),
-            postsmoother = AlgebraicMultigrid.Jacobi(rand(size(W,
-                1)))))
+            presmoother = AlgebraicMultigrid.Jacobi(rand(size(W, 1))),
+            postsmoother = AlgebraicMultigrid.Jacobi(rand(size(W, 1)))))
     end
 end
 
diff --git a/docs/src/tutorials/bvp_example.md b/docs/src/tutorials/bvp_example.md
index 5fb93439a..5402b9733 100644
--- a/docs/src/tutorials/bvp_example.md
+++ b/docs/src/tutorials/bvp_example.md
@@ -15,7 +15,7 @@ g(u) &= \vec{0}
 
 ## Example 1: Simple Pendulum
 
-The concrete example that we are solving is the simple pendulum ``\ddot{u}+\frac{g}{L}sin(u)=0`` on the time interval ``t\in[0,\frac{\pi}{2}]``. First, we need to define the ODE
+The concrete example that we are solving is the simple pendulum ``\ddot{u}+\frac{g}{L}\sin(u)=0`` on the time interval ``t\in[0,\frac{\pi}{2}]``. First, we need to define the ODE
 
 ```@example bvp
 import BoundaryValueDiffEq as BVP
@@ -36,7 +36,7 @@ There are two problem types available:
   - A problem type for general boundary conditions `BVProblem` (including conditions that may be anywhere/ everywhere on the integration interval, aka multi-points BVP).
   - A problem type for boundaries that are specified at the beginning and the end of the integration interval `TwoPointBVProblem`(aka two-points BVP)
 
-The boundary conditions are specified by a function that calculates the residual in-place from the problem solution, such that the residual is $\vec{0}$ when the boundary condition is satisfied.
+The boundary conditions are specified by a function that calculates the residual in-place from the problem solution, such that the residual is ``\vec{0}`` when the boundary condition is satisfied.
 
 There are collocation and shooting methods for addressing boundary value problems in DifferentialEquations.jl. We need to use appropriate [available BVP solvers](@ref bvp_solve) to solve `BVProblem`. In this example, we use `MIRK4` to solve the simple pendulum example.
 
@@ -100,26 +100,29 @@ Suppose we want to solve the second order BVP system which can be formulated as
 ```math
 \begin{cases}
 u_1''(x)= u_2(x),\\
-\epsilon u_2''(x)=-u_1(x)u_2'(x)- u_3(x)u_3'(x),\\
-\epsilon u_3''(x)=u_1'(x)u_3(x)- u_1(x) u_3 '(x)
+ε u_2''(x)=-u_1(x)u_2'(x)- u_3(x)u_3'(x),\\
+ε u_3''(x)=u_1'(x)u_3(x)- u_1(x) u_3 '(x)
 \end{cases}
 ```
 
 with initial conditions:
 
 ```math
-u_1(0) = u_1'(0)= u_1(1)=u_1'(1)=0,u_3(0)=
--1, u_3(1)=1
+\begin{align*}
+u_1(0) &= u_1'(0)= u_1(1)=u_1'(1)=0, \\
+u_3(0) &= -1, \\
+u_3(1) &= 1
+\end{align*}
 ```
 
 The common way of solving the second order BVP is to define intermediate variables and transform the second order system into first order one, however, DifferentialEquations.jl allows the direct solving of second order BVP system to achieve more efficiency and higher continuity of the numerical solution.
 
 ```@example bvp
 function f!(ddu, du, u, p, t)
-    ϵ = 0.1
+    ε = 0.1
     ddu[1] = u[2]
-    ddu[2] = (-u[1] * du[2] - u[3] * du[3]) / ϵ
-    ddu[3] = (du[1] * u[3] - u[1] * du[3]) / ϵ
+    ddu[2] = (-u[1] * du[2] - u[3] * du[3]) / ε
+    ddu[3] = (du[1] * u[3] - u[1] * du[3]) / ε
 end
 function bc!(res, du, u, p, t)
     res[1] = u(0.0)[1]
@@ -142,7 +145,7 @@ Consider a semi-explicit boundary value differential-algebraic equation formulat
 
 ```math
 \begin{cases}
-x_1'=(\epsilon+x_2-p_2(t))y+p_1'(t) \\
+x_1'=(ε+x_2-p_2(t))y+p_1'(t) \\
 x_2'=p_2'(t) \\
 x_3'=y \\
 0=(x_1-p_1(t))(y-e^t)
@@ -152,18 +155,21 @@ x_3'=y \\
 with boundary conditions
 
 ```math
-x_1(0)=0,x_3(0)=1,x_2(1)=\sin(1)
+\begin{align*}
+x_1(0)&=0, \\
+x_3(0)&=1, \\
+x_2(1)&=\sin(1)
+\end{align*}
 ```
 
 We need to choose the Ascher methods for solving BVDAEs.
 
 ```@example bvp
 function f!(du, u, p, t)
-    e = 2.7
     du[1] = (1 + u[2] - sin(t)) * u[4] + cos(t)
     du[2] = cos(t)
     du[3] = u[4]
-    du[4] = (u[1] - sin(t)) * (u[4] - e^t)
+    du[4] = (u[1] - sin(t)) * (u[4] - exp(t))
 end
 function bc!(res, u, p, t)
     res[1] = u[1]
diff --git a/docs/src/tutorials/dae_example.md b/docs/src/tutorials/dae_example.md
index 1b236ba4e..bf909853c 100644
--- a/docs/src/tutorials/dae_example.md
+++ b/docs/src/tutorials/dae_example.md
@@ -25,8 +25,8 @@ In previous tutorials, we wrote this equation as:
 ```math
 \begin{aligned}
 dy_1 &= -0.04 y_1 + 10^4 y_2 y_3 \\
-dy_2 &=  0.04 y_1 - 10^4 y_2 y_3 - 3*10^7 y_{2}^2 \\
-dy_3 &= 3*10^7 y_{2}^2 \\
+dy_2 &=  0.04 y_1 - 10^4 y_2 y_3 - 3×10^7 y_{2}^2 \\
+dy_3 &= 3×10^7 y_{2}^2 \\
 \end{aligned}
 ```
 
@@ -35,7 +35,7 @@ But we can instead write this with a conservation relation:
 ```math
 \begin{aligned}
 \frac{dy_1}{dt} &= -0.04 y_1 + 10^4 y_2 y_3 \\
-\frac{dy_2}{dt} &=  0.04 y_1 - 10^4 y_2 y_3 - 3*10^7 y_{2}^2 \\
+\frac{dy_2}{dt} &=  0.04 y_1 - 10^4 y_2 y_3 - 3×10^7 y_{2}^2 \\
 1 &=  y_{1} + y_{2} + y_{3} \\
 \end{aligned}
 ```
@@ -94,7 +94,7 @@ equation. The Robertson model can be written in the form:
 ```math
 \begin{aligned}
 \frac{dy_1}{dt} &= -0.04y₁ + 10^4 y_2 y_3 \\
-\frac{dy_2}{dt} &= 0.04 y_1 - 10^4 y_2 y_3 - 3*10^7 y_{2}^2 \\
+\frac{dy_2}{dt} &= 0.04 y_1 - 10^4 y_2 y_3 - 3×10^7 y_{2}^2 \\
 1 &=  y_{1} + y_{2} + y_{3} \\
 \end{aligned}
 ```
diff --git a/docs/src/tutorials/faster_ode_example.md b/docs/src/tutorials/faster_ode_example.md
index da0c669b7..ba9281836 100644
--- a/docs/src/tutorials/faster_ode_example.md
+++ b/docs/src/tutorials/faster_ode_example.md
@@ -216,8 +216,8 @@ ROBER):
 ```math
 \begin{aligned}
 \frac{dy_1}{dt} &= -0.04y₁ + 10^4 y_2 y_3 \\
-\frac{dy_2}{dt} &= 0.04 y_1 - 10^4 y_2 y_3 - 3*10^7 y_{2}^2 \\
-\frac{dy_3}{dt} &= 3*10^7 y_{2}^2 \\
+\frac{dy_2}{dt} &= 0.04 y_1 - 10^4 y_2 y_3 - 3×10^7 y_{2}^2 \\
+\frac{dy_3}{dt} &= 3×10^7 y_{2}^2 \\
 \end{aligned}
 ```
 
@@ -380,10 +380,10 @@ Let's optimize the solution of a Reaction-Diffusion PDE's discretization.
 In its discretized form, this is the ODE:
 
 ```math
-\begin{align}
+\begin{align*}
 du &= D_1 (A_y u + u A_x) + \frac{au^2}{v} + \bar{u} - \alpha u\\
 dv &= D_2 (A_y v + v A_x) + a u^2 + \beta v
-\end{align}
+\end{align*}
 ```
 
 where ``u``, ``v``, and ``A`` are matrices. Here, we will use the simplified
@@ -684,9 +684,9 @@ The last thing to do is then ***optimize our algorithm choice***. We have been
 using `DE.Tsit5()` as our test algorithm, but in reality this problem is a stiff
 PDE discretization and thus one recommendation is to use `Sundials.CVODE_BDF()`. However,
 instead of using the default dense Jacobian, we should make use of the sparse
-Jacobian afforded by the problem. The Jacobian is the matrix $\frac{df_i}{dr_j}$,
-where $r$ is read by the linear index (i.e. down columns). But since the $u$
-variables depend on the $v$, the band size here is large, and thus this will
+Jacobian afforded by the problem. The Jacobian is the matrix ``\frac{df_i}{dr_j}``,
+where ``r`` is read by the linear index (i.e. down columns). But since the ``u``
+variables depend on the ``v``, the band size here is large, and thus this will
 not do well with a Banded Jacobian solver. Instead, we utilize sparse Jacobian
 algorithms. `Sundials.CVODE_BDF` allows us to use a sparse Newton-Krylov solver by
 setting `linear_solver = :GMRES`.
diff --git a/docs/src/tutorials/sde_example.md b/docs/src/tutorials/sde_example.md
index e984f8e92..fef3d1ce2 100644
--- a/docs/src/tutorials/sde_example.md
+++ b/docs/src/tutorials/sde_example.md
@@ -28,7 +28,7 @@ import StochasticDiffEq as SDE
 u₀ = 1 / 2
 f(u, p, t) = α * u
 g(u, p, t) = β * u
-dt = 1 // 2^(4)
+dt = 1 // 2^4
 tspan = (0.0, 1.0)
 prob = SDE.SDEProblem(f, g, u₀, tspan)
 ```
@@ -132,7 +132,7 @@ nothing # hide
 
 !!! warning
     
-    If you use a custom noise process, you might need to specify it in a custom prob_func
+    If you use a custom noise process, you might need to specify it in a custom `prob_func`
     in the EnsembleProblem constructor, as each trajectory needs its own noise process.
 
 Many more controls are defined at the [Ensemble simulations page](@ref ensemble),
@@ -311,7 +311,7 @@ prob = SDE.SDEProblem(f!, g!, ones(2), (0.0, 1.0), noise_rate_prototype = A)
 and now `g!(u,p,t)` writes into a sparse matrix, and `g!(u,p,t)*dW` is sparse matrix
 multiplication.
 
-## Example 4: Colored Noise
+## Example 5: Colored Noise
 
 Colored noise can be defined [using the Noise Process interface](@ref noise_process).
 In that portion of the docs, it is shown how to define your own noise process
@@ -329,9 +329,11 @@ This is discussed [in the SDE solvers page](@ref sde_solve).
 Let's define the Heston equation from financial mathematics:
 
 ```math
-dS = μSdt + \sqrt{v}SdW_1 \\
-dv = κ(Θ-v)dt + σ\sqrt{v}dW_2 \\
-dW_1 dW_2 = ρ dt
+\begin{align*}
+dS &= μS \, dt + \sqrt{v}S \, dW_1 \\
+dv &= κ(Θ-v) \, dt + σ\sqrt{v} \, dW_2 \\
+dW_1 \, dW_2 &= ρ \, dt
+\end{align*}
 ```
 
 In this problem, we have a diagonal noise problem given by: