diff --git a/tutorials-v4/lectures/Lecture-0-Introduction-to-QuTiP.md b/tutorials-v4/lectures/Lecture-0-Introduction-to-QuTiP.md
index a5ade5ef..3f9e9c58 100644
--- a/tutorials-v4/lectures/Lecture-0-Introduction-to-QuTiP.md
+++ b/tutorials-v4/lectures/Lecture-0-Introduction-to-QuTiP.md
@@ -16,7 +16,7 @@ jupyter:
 
 Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
 
-This lecture series was developed by J.R. Johannson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
+This lecture series was developed by J.R. Johansson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
 
 This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the [qutip-tutorials repository](https://github.com/qutip/qutip-tutorials). This lecture and other tutorial notebooks are indexed at the [QuTiP Tutorial webpage](https://qutip.org/tutorials.html).
 
diff --git a/tutorials-v4/lectures/Lecture-1-Jaynes-Cumming-model.md b/tutorials-v4/lectures/Lecture-1-Jaynes-Cumming-model.md
index 8991dec6..60d59a4a 100644
--- a/tutorials-v4/lectures/Lecture-1-Jaynes-Cumming-model.md
+++ b/tutorials-v4/lectures/Lecture-1-Jaynes-Cumming-model.md
@@ -17,7 +17,7 @@ jupyter:
 
 Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
 
-This lecture series was developed by J.R. Johannson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
+This lecture series was developed by J.R. Johansson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
 
 This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the [qutip-tutorials repository](https://github.com/qutip/qutip-tutorials). This lecture and other tutorial notebooks are indexed at the [QuTiP Tutorial webpage](https://qutip.org/tutorials.html).
 
@@ -65,7 +65,7 @@ tlist = np.linspace(0, 25, 101)
 ### Setup the operators, the Hamiltonian and initial state
 
 ```python
-# intial state
+# initial state
 psi0 = tensor(basis(N, 0), basis(2, 1))  # start with an excited atom
 
 # operators
@@ -132,7 +132,7 @@ axes.set_title("Vacuum Rabi oscillations");
 
 In addition to the cavity's and atom's excitation probabilities, we may also be interested in for example the wigner function as a function of time. The Wigner function can give some valuable insight in the nature of the state of the resonators. 
 
-To calculate the Wigner function in QuTiP, we first recalculte the evolution without specifying any expectation value operators, which will result in that the solver return a list of density matrices for the system for the given time coordinates.
+To calculate the Wigner function in QuTiP, we first recalculate the evolution without specifying any expectation value operators, which will result in that the solver returns a list of density matrices for the system for the given time coordinates.
 
 ```python
 output = mesolve(H, psi0, tlist, c_ops, [])
@@ -163,7 +163,7 @@ For each of these points in time we need to:
 
  1. Find the system density matrix for the points in time that we are interested in.
  2. Trace out the atom and obtain the reduced density matrix for the cavity.
- 3. Calculate and visualize the Wigner function fo the reduced cavity density matrix.
+ 3. Calculate and visualize the Wigner function for the reduced cavity density matrix.
 
 ```python
 # find the indices of the density matrices for the times we are interested in
@@ -206,7 +206,7 @@ for idx, rho in enumerate(rho_list):
     axes[idx].set_title(r"$t = %.1f$" % tlist[t_idx][idx], fontsize=16)
 ```
 
-At $t =0$, the cavity is in it's ground state. At $t = 5, 15, 25$ it reaches it's maxium occupation in this Rabi-vacuum oscillation process. We can note that for $t=5$ and $t=15$ the Wigner function has negative values, indicating a truely quantum mechanical state. At $t=25$, however, the wigner function no longer has negative values and can therefore be considered a classical state.
+At $t =0$, the cavity is in its ground state. At $t = 5, 15, 25$ it reaches its maximum occupation in this Rabi-vacuum oscillation process. We can note that for $t=5$ and $t=15$ the Wigner function has negative values, indicating a truly quantum mechanical state. At $t=25$, however, the Wigner function no longer has negative values and can therefore be considered a classical state.
 
 
 ### Alternative view of the same thing
diff --git a/tutorials-v4/lectures/Lecture-10-cQED-dispersive-regime.md b/tutorials-v4/lectures/Lecture-10-cQED-dispersive-regime.md
index cdffe32d..ae3ee72a 100644
--- a/tutorials-v4/lectures/Lecture-10-cQED-dispersive-regime.md
+++ b/tutorials-v4/lectures/Lecture-10-cQED-dispersive-regime.md
@@ -16,7 +16,7 @@ jupyter:
 
 Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
 
-This lecture series was developed by J.R. Johannson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
+This lecture series was developed by J.R. Johansson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
 
 This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the [qutip-tutorials repository](https://github.com/qutip/qutip-tutorials). This lecture and other tutorial notebooks are indexed at the [QuTiP Tutorial webpage](https://qutip.org/tutorials.html).
 
diff --git a/tutorials-v4/lectures/Lecture-11-Charge-Qubits.md b/tutorials-v4/lectures/Lecture-11-Charge-Qubits.md
index 70ed1617..c3aa172c 100644
--- a/tutorials-v4/lectures/Lecture-11-Charge-Qubits.md
+++ b/tutorials-v4/lectures/Lecture-11-Charge-Qubits.md
@@ -16,7 +16,7 @@ jupyter:
 
 Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
 
-This lecture series was developed by J.R. Johannson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
+This lecture series was developed by J.R. Johansson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
 
 This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the [qutip-tutorials repository](https://github.com/qutip/qutip-tutorials). This lecture and other tutorial notebooks are indexed at the [QuTiP Tutorial webpage](https://qutip.org/tutorials.html).
 
@@ -39,7 +39,7 @@ where $E_C$ is the charge energy, $E_J$ is the Josephson energy, and $\left| n\r
 
 #### References
 
- * [J. Koch et al, Phys. Rec. A 76, 042319 (2007)](https://journals.aps.org/pra/abstract/10.1103/PhysRevA.76.042319)
+ * [J. Koch et al, Phys. Rev. A 76, 042319 (2007)](https://journals.aps.org/pra/abstract/10.1103/PhysRevA.76.042319)
  * [Y.A. Pashkin et al, Quantum Inf Process 8, 55 (2009)](http://dx.doi.org/10.1007/s11128-009-0101-5)
 
 
@@ -176,7 +176,7 @@ energies = np.array([hamiltonian(Ec, Ej, N, ng).eigenenergies()
 plot_energies(ng_vec, energies, ymax=(50, 3));
 ```
 
-Note that the energy-level splitting is essentially independent of the gate bias $n_g$, at least for the lowest few states. This device insensitive to charge noise. But at the same time the two lowest energy states are no longer well separated from higher states (it has become more like an harmonic oscillator). But some anharmonicity still remains, and it can still be used as a qubit if the leakage of occupation probability of the higher states can be kept under control.
+Note that the energy-level splitting is essentially independent of the gate bias $n_g$, at least for the lowest few states. This device is insensitive to charge noise. But at the same time the two lowest energy states are no longer well separated from higher states (it has become more like a harmonic oscillator). But some anharmonicity still remains, and it can still be used as a qubit if the leakage of occupation probability to the higher states can be kept under control.
 
 
 ## Focus on the two lowest energy states
diff --git a/tutorials-v4/lectures/Lecture-12-Decay-into-a-squeezed-vacuum-field.md b/tutorials-v4/lectures/Lecture-12-Decay-into-a-squeezed-vacuum-field.md
index 31168ea3..7f1f1af7 100644
--- a/tutorials-v4/lectures/Lecture-12-Decay-into-a-squeezed-vacuum-field.md
+++ b/tutorials-v4/lectures/Lecture-12-Decay-into-a-squeezed-vacuum-field.md
@@ -17,7 +17,7 @@ jupyter:
 
 Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
 
-This lecture series was developed by J.R. Johannson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
+This lecture series was developed by J.R. Johansson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
 
 This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the [qutip-tutorials repository](https://github.com/qutip/qutip-tutorials). This lecture and other tutorial notebooks are indexed at the [QuTiP Tutorial webpage](https://qutip.org/tutorials.html).
 
@@ -70,7 +70,7 @@ w_th = 0.0 * 2 * np.pi
 
 ```python
 # the number of average excitations in the
-# environment mode w0 at temperture w_th
+# environment mode w0 at temperature w_th
 Nth = n_thermal(w0, w_th)
 
 Nth
@@ -126,7 +126,7 @@ L0 = liouvillian(H, c_ops)
 L0
 ```
 
-Next we manually construct the Liouvillian for the effect of the squeeing in the environment, which is not on standard form we can therefore not use the `liouvillian` function in QuTiP
+Next we manually construct the Liouvillian for the effect of the squeezing in the environment, which is not in standard form so we can therefore not use the `liouvillian` function in QuTiP
 
 ```python
 Lsq = -gamma0 * M * spre(sp) * spost(sp) - gamma0 * \
diff --git a/tutorials-v4/lectures/Lecture-13-Resonance-flourescence.md b/tutorials-v4/lectures/Lecture-13-Resonance-flourescence.md
index fc269284..f74f29bf 100644
--- a/tutorials-v4/lectures/Lecture-13-Resonance-flourescence.md
+++ b/tutorials-v4/lectures/Lecture-13-Resonance-flourescence.md
@@ -12,12 +12,12 @@ jupyter:
     name: python3
 ---
 
-# Lecture 13 - Resonance flourescence
+# Lecture 13 - Resonance fluorescence
 
 
 Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
 
-This lecture series was developed by J.R. Johannson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
+This lecture series was developed by J.R. Johansson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
 
 This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the [qutip-tutorials repository](https://github.com/qutip/qutip-tutorials). This lecture and other tutorial notebooks are indexed at the [QuTiP Tutorial webpage](https://qutip.org/tutorials.html).
 
@@ -34,6 +34,10 @@ from qutip import (about, basis, correlation_2op_1t, mesolve, n_thermal, num,
 <!-- #region -->
 ## Introduction
 
+Resonance fluorescence, is the interaction between a two-level quantum system and a coherent driving field.
+
+The Hamiltonian for a two-level system under resonant driving, is given by:
+
 $\displaystyle H_L = -\frac{\Omega}{2}(\sigma_+ + \sigma_-)$
 
 
@@ -44,12 +48,18 @@ $\displaystyle \frac{d}{dt}\rho = -i[H_L, \rho] + \gamma_0(N+1)\left(\sigma_-\rh
 ### Problem definition in QuTiP
 
 ```python
+# Rabi frequency (strength of the driving field)
 Omega = 1.0 * 2 * np.pi
 ```
 
 ```python
+# decay rate of the two-level system
 gamma0 = 0.05
+
+# thermal frequency of the bath
 w_th = 0.0
+
+# thermal occupation number for the bath
 N = n_thermal(Omega, w_th)
 ```
 
@@ -90,34 +100,41 @@ axes[0].plot(result.times, result.expect[2], "b",
 axes[0].legend()
 axes[0].set_ylim(-1, 1)
 
+axes[0].set_ylabel(r"$\langle\sigma_{x,y,z}\rangle$", fontsize=16)
+axes[0].set_title("Bloch Vector Components vs Time")
 
 axes[1].plot(result.times, result.expect[5], "b", label=r"$P_e$")
 
-# axes[1].set_ylabel(r'$\langle\sigma_z\rangle$', fontsize=16)
+axes[1].set_ylabel(r"$P_g$", fontsize=16)
 axes[1].set_xlabel("time", fontsize=16)
 axes[1].legend()
-axes[1].set_ylim(0, 1);
+axes[1].set_ylim(0, 1)
+axes[1].set_title("Ground State Population vs Time");
 ```
 
 ```python
 fig, ax = plt.subplots(1, 1, figsize=(12, 6), sharex=True)
 
-
 for idx, gamma0 in enumerate([0.1 * Omega, 0.5 * Omega, 1.0 * Omega]):
 
     HL, c_ops = system_spec(Omega, gamma0, N)
     result = mesolve(HL, psi0, tlist, c_ops, e_ops)
 
-    ax.plot(result.times, result.expect[5], "b",
-            label=r"$\langle\sigma_z\rangle$")
+    ax.plot(
+        result.times, result.expect[5], "b",
+        label=fr"$P_g$ ($\gamma_0={gamma0:.2f}$)"
+    )
 
-ax.set_ylim(0, 1);
+ax.set_ylim(0, 1)
+ax.set_xlabel("time", fontsize=16)
+ax.set_ylabel(r"$P_g$", fontsize=16)
+ax.set_title("Ground State Population for Different Decay Rates")
+ax.legend();
 ```
 
 ```python
 fig, ax = plt.subplots(1, 1, figsize=(12, 6), sharex=True)
 
-
 for idx, gamma0 in enumerate([0.1 * Omega, 0.5 * Omega, 1.0 * Omega]):
 
     HL, c_ops = system_spec(Omega, gamma0, N)
@@ -125,10 +142,14 @@ for idx, gamma0 in enumerate([0.1 * Omega, 0.5 * Omega, 1.0 * Omega]):
 
     ax.plot(
         result.times, np.imag(result.expect[4]),
-        label=r"im $\langle\sigma_+\rangle$"
+        label=fr"Im $\langle\sigma_+\rangle$ ($\gamma_0={gamma0:.2f}$)"
     )
 
-ax.set_ylim(-0.5, 0.5);
+ax.set_ylim(-0.5, 0.5)
+ax.set_xlabel("time", fontsize=16)
+ax.set_ylabel(r"Im $\langle\sigma_+\rangle$", fontsize=16)
+ax.set_title(r"Imaginary Part of $\langle\sigma_+\rangle$ for Different Decay Rates")
+ax.legend();
 ```
 
 ```python
@@ -142,12 +163,21 @@ for idx, gamma0 in enumerate([2 * Omega, 0.5 * Omega, 0.25 * Omega]):
     corr_vec = correlation_2op_1t(HL, None, taulist, c_ops, sigmap(), sigmam())
     w, S = spectrum_correlation_fft(taulist, corr_vec)
 
-    axes[0].plot(taulist, corr_vec, label=r"$<\sigma_+(\tau)\sigma_-(0)>$")
-    axes[1].plot(-w / (gamma0), S, "b", label=r"$S(\omega)$")
-    axes[1].plot(w / (gamma0), S, "b", label=r"$S(\omega)$")
+    axes[0].plot(taulist, corr_vec, label=fr"$\gamma_0={gamma0:.2f}$")
+    axes[1].plot(-w / (gamma0), S, "b", label=fr"$\gamma_0={gamma0:.2f}$")
+    axes[1].plot(w / (gamma0), S, "b")
 
 axes[0].set_xlim(0, 10)
-axes[1].set_xlim(-5, 5);
+axes[0].set_xlabel(r"$\tau$", fontsize=16)
+axes[0].set_ylabel(r"$\langle\sigma_+(\tau)\sigma_-(0)\rangle$", fontsize=16)
+axes[0].set_title("Two-Time Correlation Function")
+axes[0].legend()
+
+axes[1].set_xlim(-5, 5)
+axes[1].set_xlabel(r"$\omega/\gamma_0$", fontsize=16)
+axes[1].set_ylabel(r"$S(\omega)$", fontsize=16)
+axes[1].set_title("Resonance Fluorescence Spectrum")
+axes[1].legend();
 ```
 
 ### Software versions
diff --git a/tutorials-v4/lectures/Lecture-14-Kerr-nonlinearities.md b/tutorials-v4/lectures/Lecture-14-Kerr-nonlinearities.md
index b6bb4288..4697123b 100644
--- a/tutorials-v4/lectures/Lecture-14-Kerr-nonlinearities.md
+++ b/tutorials-v4/lectures/Lecture-14-Kerr-nonlinearities.md
@@ -16,7 +16,7 @@ jupyter:
 
 Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
 
-This lecture series was developed by J.R. Johannson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
+This lecture series was developed by J.R. Johansson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
 
 This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the [qutip-tutorials repository](https://github.com/qutip/qutip-tutorials). This lecture and other tutorial notebooks are indexed at the [QuTiP Tutorial webpage](https://qutip.org/tutorials.html).
 
@@ -40,7 +40,7 @@ The Kerr effect describes a self-interaction electromagnetic quantum field which
 
 $\displaystyle H = \frac{1}{2}\chi (a^\dagger)^2a^2$
 
-where $\chi$ is related to the third-order nonlinear suseptibility. The Kerr effect is one of the typical nonlinearities that can occur in quantum optics due to a nonlinear medium.
+where $\chi$ is related to the third-order nonlinear susceptibility. The Kerr effect is one of the typical nonlinearities that can occur in quantum optics due to a nonlinear medium.
 
 In this notebook we'll see how to setup the model in QuTiP and look at some interesting properties of the states that evolve according to this Hamiltonian.
 
@@ -182,12 +182,12 @@ psi0 = coherent(N, 2.0)
 
 ```python
 # and evolve the state under the influence of the hamiltonian.
-# by passing an empty list as expecation value operators argument,
+# by passing an empty list as expectation value operators argument,
 # we get the full state of the system in result.states
 result = mesolve(H, psi0, tlist, [], [])
 ```
 
-First, let's look at how the expecation values and variances of the photon number operator $n$ and the $x$ and $p$ quadratures evolve in time:
+First, let's look at how the expectation values and variances of the photon number operator $n$ and the $x$ and $p$ quadratures evolve in time:
 
 ```python
 plot_expect_with_variance(N, [n, x, p], [r"n", r"x", r"p"], result.states);
@@ -201,7 +201,7 @@ To verify that the photon distribution indeed is time-independent, we can plot t
 plot_fock_distribution_vs_time(tlist, result.states);
 ```
 
-So the fock state distribution is constant, but let's see how the Wigner function of the state evolves in time. To best illustrate the dynamics of the Winger function we make a short movie that show the Wigner function from time $t=0$ to the the final time of the evolution.
+So the fock state distribution is constant, but let's see how the Wigner function of the state evolves in time. To best illustrate the dynamics of the Wigner function we make a short movie that shows the Wigner function from time $t=0$ to the final time of the evolution.
 
 ```python
 fig, ax = plt.subplots(1, 1, figsize=(8, 8))
@@ -226,7 +226,7 @@ display_embedded_video("animation-kerr-coherent-state.mp4")
 
 Isn't that interesting! The dynamics is periodic, and we evolved the state for exactly one period, so that the final state is equal to the initial state.
 
-In between there is interesting stuff going on. For example, after half the period the state ends up in something that look very much like a cat-state superposition of coherent states!
+In between there is interesting stuff going on. For example, after half the period the state ends up in something that looks very much like a cat-state superposition of coherent states!
 
 ```python
 fig, ax = plt.subplots(1, 1, figsize=(8, 8))
diff --git a/tutorials-v4/lectures/Lecture-15-Nonclassically-driven-atoms.md b/tutorials-v4/lectures/Lecture-15-Nonclassically-driven-atoms.md
index 00a6e94c..e25c5e5a 100644
--- a/tutorials-v4/lectures/Lecture-15-Nonclassically-driven-atoms.md
+++ b/tutorials-v4/lectures/Lecture-15-Nonclassically-driven-atoms.md
@@ -17,7 +17,7 @@ jupyter:
 
 Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
 
-This lecture series was developed by J.R. Johannson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
+This lecture series was developed by J.R. Johansson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
 
 This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the [qutip-tutorials repository](https://github.com/qutip/qutip-tutorials). This lecture and other tutorial notebooks are indexed at the [QuTiP Tutorial webpage](https://qutip.org/tutorials.html).
 
diff --git a/tutorials-v4/lectures/Lecture-16-Gallery-of-Wigner-functions.md b/tutorials-v4/lectures/Lecture-16-Gallery-of-Wigner-functions.md
index bbfbf4cc..66346fe1 100644
--- a/tutorials-v4/lectures/Lecture-16-Gallery-of-Wigner-functions.md
+++ b/tutorials-v4/lectures/Lecture-16-Gallery-of-Wigner-functions.md
@@ -17,7 +17,7 @@ jupyter:
 
 Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
 
-This lecture series was developed by J.R. Johannson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
+This lecture series was developed by J.R. Johansson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
 
 This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the [qutip-tutorials repository](https://github.com/qutip/qutip-tutorials). This lecture and other tutorial notebooks are indexed at the [QuTiP Tutorial webpage](https://qutip.org/tutorials.html).
 
diff --git a/tutorials-v4/lectures/Lecture-2A-Cavity-Qubit-Gates.md b/tutorials-v4/lectures/Lecture-2A-Cavity-Qubit-Gates.md
index 868bb6b9..e3553cfa 100644
--- a/tutorials-v4/lectures/Lecture-2A-Cavity-Qubit-Gates.md
+++ b/tutorials-v4/lectures/Lecture-2A-Cavity-Qubit-Gates.md
@@ -16,7 +16,7 @@ jupyter:
 
 Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
 
-This lecture series was developed by J.R. Johannson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
+This lecture series was developed by J.R. Johansson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
 
 This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the [qutip-tutorials repository](https://github.com/qutip/qutip-tutorials). This lecture and other tutorial notebooks are indexed at the [QuTiP Tutorial webpage](https://qutip.org/tutorials.html).
 
@@ -68,7 +68,7 @@ sm1 = tensor(qeye(N), destroy(2), qeye(2))
 sz1 = tensor(qeye(N), sigmaz(), qeye(2))
 n1 = sm1.dag() * sm1
 
-# oeprators for qubit 2
+# operators for qubit 2
 sm2 = tensor(qeye(N), qeye(2), destroy(2))
 sz2 = tensor(qeye(N), qeye(2), sigmaz())
 n2 = sm2.dag() * sm2
diff --git a/tutorials-v4/lectures/Lecture-2B-Single-Atom-Lasing.md b/tutorials-v4/lectures/Lecture-2B-Single-Atom-Lasing.md
index f9efd6ee..a35ac9a0 100644
--- a/tutorials-v4/lectures/Lecture-2B-Single-Atom-Lasing.md
+++ b/tutorials-v4/lectures/Lecture-2B-Single-Atom-Lasing.md
@@ -17,7 +17,7 @@ jupyter:
 
 Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
 
-This lecture series was developed by J.R. Johannson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
+This lecture series was developed by J.R. Johansson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
 
 This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the [qutip-tutorials repository](https://github.com/qutip/qutip-tutorials). This lecture and other tutorial notebooks are indexed at the [QuTiP Tutorial webpage](https://qutip.org/tutorials.html).
 
@@ -90,7 +90,7 @@ tlist = np.linspace(0, 150, 101)
 ### Setup the operators, the Hamiltonian and initial state
 
 ```python
-# intial state
+# initial state
 psi0 = tensor(basis(N, 0), basis(2, 0))  # start without excitations
 
 # operators
@@ -237,7 +237,7 @@ References:
  * [S. Ashhab, J.R. Johansson, A.M. Zagoskin, F. Nori, New J. Phys. 11, 023030 (2009)](http://dx.doi.org/10.1088/1367-2630/11/2/023030)
 
 ```python
-def calulcate_avg_photons(N, Gamma):
+def calculate_avg_photons(N, Gamma):
 
     # collapse operators
     c_ops = []
@@ -278,7 +278,7 @@ n_avg_vec = []
 g2_vec = []
 
 for Gamma in Gamma_vec:
-    n_avg, g2 = calulcate_avg_photons(N, Gamma)
+    n_avg, g2 = calculate_avg_photons(N, Gamma)
     n_avg_vec.append(n_avg)
     g2_vec.append(g2)
 ```
diff --git a/tutorials-v4/lectures/Lecture-3A-Dicke-model.md b/tutorials-v4/lectures/Lecture-3A-Dicke-model.md
index 3aaae2ab..1049fce4 100644
--- a/tutorials-v4/lectures/Lecture-3A-Dicke-model.md
+++ b/tutorials-v4/lectures/Lecture-3A-Dicke-model.md
@@ -16,7 +16,7 @@ jupyter:
 
 Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
 
-This lecture series was developed by J.R. Johannson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
+This lecture series was developed by J.R. Johansson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
 
 This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the [qutip-tutorials repository](https://github.com/qutip/qutip-tutorials). This lecture and other tutorial notebooks are indexed at the [QuTiP Tutorial webpage](https://qutip.org/tutorials.html).
 
@@ -201,7 +201,7 @@ fig.tight_layout()
 * [Lambert et al., Phys. Rev. Lett. 92, 073602 (2004)](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.92.073602).
 
 ```python
-def calulcate_entropy(M, N, g_vec):
+def calculate_entropy(M, N, g_vec):
 
     j = N / 2.0
     n = 2 * j + 1
@@ -241,7 +241,7 @@ fig, axes = plt.subplots(1, 1, figsize=(12, 6))
 
 for NN in N_vec:
 
-    entropy_cavity, entropy_spin = calulcate_entropy(MM, NN, g_vec)
+    entropy_cavity, entropy_spin = calculate_entropy(MM, NN, g_vec)
 
     axes.plot(g_vec, entropy_cavity, "b", label="N = %d" % NN)
     axes.plot(g_vec, entropy_spin, "r--")
diff --git a/tutorials-v4/lectures/Lecture-3B-Jaynes-Cumming-model-with-ultrastrong-coupling.md b/tutorials-v4/lectures/Lecture-3B-Jaynes-Cumming-model-with-ultrastrong-coupling.md
index cce8fc10..d68aeff2 100644
--- a/tutorials-v4/lectures/Lecture-3B-Jaynes-Cumming-model-with-ultrastrong-coupling.md
+++ b/tutorials-v4/lectures/Lecture-3B-Jaynes-Cumming-model-with-ultrastrong-coupling.md
@@ -17,7 +17,7 @@ jupyter:
 
 Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
 
-This lecture series was developed by J.R. Johannson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
+This lecture series was developed by J.R. Johansson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
 
 This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the [qutip-tutorials repository](https://github.com/qutip/qutip-tutorials). This lecture and other tutorial notebooks are indexed at the [QuTiP Tutorial webpage](https://qutip.org/tutorials.html).
 
@@ -104,7 +104,7 @@ for g in g_vec:
 
     H = H0 + g * H1
 
-    # find the groundstate and its energy
+    # find the ground state and its energy
     gnd_energy, gnd_state = H.groundstate()
 
     # store the ground state
@@ -118,7 +118,7 @@ na_expt = expect(na, psi_list)  # qubit  occupation probability
 nc_expt = expect(nc, psi_list)  # cavity occupation probability
 ```
 
-Plot the ground state occupation probabilities of the cavity and the atom as a function of coupling strenght. Note that for large coupling strength (the ultrastrong coupling regime, where $g > \omega_a,\omega_c$), the ground state has both photonic and atomic excitations.
+Plot the ground state occupation probabilities of the cavity and the atom as a function of coupling strength. Note that for large coupling strength (the ultrastrong coupling regime, where $g > \omega_a,\omega_c$), the ground state has both photonic and atomic excitations.
 
 ```python
 fig, axes = plt.subplots(1, 1, sharex=True, figsize=(8, 4))
@@ -126,7 +126,7 @@ fig, axes = plt.subplots(1, 1, sharex=True, figsize=(8, 4))
 axes.plot(g_vec / (2 * np.pi), nc_expt, "r", linewidth=2, label="cavity")
 axes.plot(g_vec / (2 * np.pi), na_expt, "b", linewidth=2, label="atom")
 axes.set_ylabel("Occupation probability", fontsize=16)
-axes.set_xlabel("coupling strenght", fontsize=16)
+axes.set_xlabel("coupling strength", fontsize=16)
 axes.legend(loc=0)
 
 fig.tight_layout()
diff --git a/tutorials-v4/lectures/Lecture-4-Correlation-Functions.md b/tutorials-v4/lectures/Lecture-4-Correlation-Functions.md
index fb585463..7795ecef 100644
--- a/tutorials-v4/lectures/Lecture-4-Correlation-Functions.md
+++ b/tutorials-v4/lectures/Lecture-4-Correlation-Functions.md
@@ -16,7 +16,7 @@ jupyter:
 
 Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
 
-This lecture series was developed by J.R. Johannson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
+This lecture series was developed by J.R. Johansson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
 
 This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the [qutip-tutorials repository](https://github.com/qutip/qutip-tutorials). This lecture and other tutorial notebooks are indexed at the [QuTiP Tutorial webpage](https://qutip.org/tutorials.html).
 
diff --git a/tutorials-v4/lectures/Lecture-5-Parametric-Amplifier.md b/tutorials-v4/lectures/Lecture-5-Parametric-Amplifier.md
index 134a6b39..dc3b7d39 100644
--- a/tutorials-v4/lectures/Lecture-5-Parametric-Amplifier.md
+++ b/tutorials-v4/lectures/Lecture-5-Parametric-Amplifier.md
@@ -16,7 +16,7 @@ jupyter:
 
 Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
 
-This lecture series was developed by J.R. Johannson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
+This lecture series was developed by J.R. Johansson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
 
 This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the [qutip-tutorials repository](https://github.com/qutip/qutip-tutorials). This lecture and other tutorial notebooks are indexed at the [QuTiP Tutorial webpage](https://qutip.org/tutorials.html).
 
@@ -95,7 +95,7 @@ for idx, psi in enumerate(output.states):
     nb_e[idx] = expect(nb, psi)
     nb_s[idx] = expect(nb * nb, psi)
 
-# substract the average squared to obtain variances
+# subtract the average squared to obtain variances
 na_s = na_s - na_e**2
 nb_s = nb_s - nb_e**2
 ```
diff --git a/tutorials-v4/lectures/Lecture-6-Quantum-Monte-Carlo-Trajectories.md b/tutorials-v4/lectures/Lecture-6-Quantum-Monte-Carlo-Trajectories.md
index be80faa8..51717533 100644
--- a/tutorials-v4/lectures/Lecture-6-Quantum-Monte-Carlo-Trajectories.md
+++ b/tutorials-v4/lectures/Lecture-6-Quantum-Monte-Carlo-Trajectories.md
@@ -17,7 +17,7 @@ jupyter:
 
 Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
 
-This lecture series was developed by J.R. Johannson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
+This lecture series was developed by J.R. Johansson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
 
 This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the [qutip-tutorials repository](https://github.com/qutip/qutip-tutorials). This lecture and other tutorial notebooks are indexed at the [QuTiP Tutorial webpage](https://qutip.org/tutorials.html).
 
@@ -125,13 +125,13 @@ The expectation values of $a^\dagger a$ are now available in array ``mc.expect[i
 
 ## Lindblad master-equation simulation and steady state
 
-For comparison with the averages of single quantum trajectories provided by the Monte-Carlo solver we here also calculate the dynamics of the Lindblad master equation, which should agree with the Monte-Carlo simultions for infinite number of trajectories.
+For comparison with the averages of single quantum trajectories provided by the Monte-Carlo solver we here also calculate the dynamics of the Lindblad master equation, which should agree with the Monte-Carlo simulations for infinite number of trajectories.
 
 ```python
 # run master equation to get ensemble average expectation values
 me = mesolve(H, psi0, tlist, c_op_list, [a.dag() * a])
 
-# calulate final state using steadystate solver
+# calculate final state using steadystate solver
 final_state = steadystate(H, c_op_list)  # find steady-state
 # find expectation value for particle number
 fexpt = expect(a.dag() * a, final_state)
diff --git a/tutorials-v4/lectures/Lecture-7-iSWAP-gate.md b/tutorials-v4/lectures/Lecture-7-iSWAP-gate.md
index ece8e269..e873a654 100644
--- a/tutorials-v4/lectures/Lecture-7-iSWAP-gate.md
+++ b/tutorials-v4/lectures/Lecture-7-iSWAP-gate.md
@@ -16,7 +16,7 @@ jupyter:
 
 Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
 
-This lecture series was developed by J.R. Johannson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
+This lecture series was developed by J.R. Johansson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
 
 This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the [qutip-tutorials repository](https://github.com/qutip/qutip-tutorials). This lecture and other tutorial notebooks are indexed at the [QuTiP Tutorial webpage](https://qutip.org/tutorials.html).
 
@@ -37,7 +37,7 @@ $\displaystyle H = g \left(\sigma_x\otimes\sigma_x + \sigma_y\otimes\sigma_y\rig
 
 where $g$ is the coupling strength. Under ideal conditions this coupling realizes the $i$-SWAP gate between the two qubit states. 
 
-Here we will solve for the dynamics of the two qubits subject to this Hamiltonian, and look at the deterioating effects of adding decoherence. We will use process tomography to visualize the gate.
+Here we will solve for the dynamics of the two qubits subject to this Hamiltonian, and look at the deteriorating effects of adding decoherence. We will use process tomography to visualize the gate.
 
 
 ### Parameters
diff --git a/tutorials-v4/lectures/Lecture-8-Adiabatic-quantum-computing.md b/tutorials-v4/lectures/Lecture-8-Adiabatic-quantum-computing.md
index bea1bd9a..1b963287 100644
--- a/tutorials-v4/lectures/Lecture-8-Adiabatic-quantum-computing.md
+++ b/tutorials-v4/lectures/Lecture-8-Adiabatic-quantum-computing.md
@@ -16,20 +16,20 @@ jupyter:
 
 Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
 
-This lecture series was developed by J.R. Johannson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
+This lecture series was developed by J.R. Johansson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
 
 This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the [qutip-tutorials repository](https://github.com/qutip/qutip-tutorials). This lecture and other tutorial notebooks are indexed at the [QuTiP Tutorial webpage](https://qutip.org/tutorials.html).
 
 
 ## Introduction
 
-In adiabatic quantum computing, an easy to prepare ground state of a Hamiltonian $H_0$ is prepared, and then the Hamiltonian is gradually transformed into $H_1$, which is constructed in such a way that the groundstate of $H_1$ encodes the solution to a difficult problem. The transformation of $H_0$ to $H_1$ can for example be written on the form
+In adiabatic quantum computing, an easy to prepare ground state of a Hamiltonian $H_0$ is prepared, and then the Hamiltonian is gradually transformed into $H_1$, which is constructed in such a way that the ground state of $H_1$ encodes the solution to a difficult problem. The transformation of $H_0$ to $H_1$ can for example be written on the form
 
 $\displaystyle H(t) = \lambda(t) H_0 + (1 - \lambda(t)) H_1$
 
-where $\lambda(t)$ is a function that goes from goes from $0$ to $1$ when $t$ goes from $0$ to $t_{\rm final}$.
+where $\lambda(t)$ is a function that goes from $0$ to $1$ when $t$ goes from $0$ to $t_{\rm final}$.
 
-If this gradual tranformation is slow enough (satisfying the adiabicity critera), the evolution of the system will remain in its ground state.
+If this gradual transformation is slow enough (satisfying the adiabaticity criteria), the evolution of the system will remain in its ground state.
 
 If the Hamiltonian is transformed from $H_0$ to $H_1$ too quickly, the system will get excited from the ground state the adiabatic computing algorithm fails.
 
@@ -194,7 +194,7 @@ for idx in range(len(taulist) - 1):
 axes[0].set_xlabel(r"$\tau$")
 axes[0].set_ylabel("Eigenenergies")
 axes[0].set_title(
-    "Energyspectrum (%d lowest values) of a chain of %d spins.\n " % (M, N)
+    "Energy spectrum (%d lowest values) of a chain of %d spins.\n " % (M, N)
     + "The occupation probabilities are encoded in the red line widths."
 )
 
diff --git a/tutorials-v4/lectures/Lecture-9-Squeezed-states-of-harmonic-oscillator.md b/tutorials-v4/lectures/Lecture-9-Squeezed-states-of-harmonic-oscillator.md
index 53f4b817..77167f96 100644
--- a/tutorials-v4/lectures/Lecture-9-Squeezed-states-of-harmonic-oscillator.md
+++ b/tutorials-v4/lectures/Lecture-9-Squeezed-states-of-harmonic-oscillator.md
@@ -16,7 +16,7 @@ jupyter:
 
 Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
 
-This lecture series was developed by J.R. Johannson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
+This lecture series was developed by J.R. Johansson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
 
 This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the [qutip-tutorials repository](https://github.com/qutip/qutip-tutorials). This lecture and other tutorial notebooks are indexed at the [QuTiP Tutorial webpage](https://qutip.org/tutorials.html).
 
@@ -35,9 +35,9 @@ from qutip import (about, basis, coherent, destroy, displace, expect, mesolve,
 
 ## Introduction
 
-In quantum mechanics, each measurement of an observable (corresponding to an Hermitian operator) results in stochastic outcome that follows some probability distribution. The expectation value of the operator is the average of many measurement outcomes, and the standard deviation of the operator describes the uncertainty in the outcomes.
+In quantum mechanics, each measurement of an observable (corresponding to a Hermitian operator) results in stochastic outcome that follows some probability distribution. The expectation value of the operator is the average of many measurement outcomes, and the standard deviation of the operator describes the uncertainty in the outcomes.
 
-This uncertainty is intrinsic in quantum mechanics, and cannot be eliminated. The Heisenberg uncertainty principle describes the minumum uncertainly for pairs of noncommuting operators. For example, the operators such $x$ and $p$, which satisfy the commutation relation $[x, p] = i\hbar$, must always satisfy $(\Delta x) (\Delta p) >= \hbar/2$ .
+This uncertainty is intrinsic in quantum mechanics, and cannot be eliminated. The Heisenberg uncertainty principle describes the minimum uncertainty for pairs of noncommuting operators. For example, the operators such $x$ and $p$, which satisfy the commutation relation $[x, p] = i\hbar$, must always satisfy $(\Delta x) (\Delta p) >= \hbar/2$ .
 
 A state that satisfies
 
@@ -47,7 +47,7 @@ is called a minimum uncertainty state, and a state for which, for example,
 
 $(\Delta x)^2 < \hbar/2$ 
 
-is called a squeezed state. Note that in this case $(\Delta p)^2$ must be larger than $\hbar/2(\Delta x)^2$ for the Heisenberg relation to hold. Squeezing a quantum state so that the variance of one operator $x$ is reduced below the minimum uncertainty limit therefore necessarily amplify the variance of operators that do not commute with $x$, such as $p$.
+is called a squeezed state. Note that in this case $(\Delta p)^2$ must be larger than $\hbar/2(\Delta x)^2$ for the Heisenberg relation to hold. Squeezing a quantum state so that the variance of one operator $x$ is reduced below the minimum uncertainty limit therefore necessarily amplifies the variance of operators that do not commute with $x$, such as $p$.
 
 For harmonic modes, squeezing of $x$ or $p$ is called quadrature squeezing, and it is probably the most common form of squeezing. 
 
diff --git a/tutorials-v5/lectures/Lecture-0-Introduction-to-QuTiP.md b/tutorials-v5/lectures/Lecture-0-Introduction-to-QuTiP.md
index f1c4a28c..49ad841d 100644
--- a/tutorials-v5/lectures/Lecture-0-Introduction-to-QuTiP.md
+++ b/tutorials-v5/lectures/Lecture-0-Introduction-to-QuTiP.md
@@ -16,7 +16,7 @@ jupyter:
 
 Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
 
-This lecture series was developed by J.R. Johannson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
+This lecture series was developed by J.R. Johansson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
 
 This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the [qutip-tutorials repository](https://github.com/qutip/qutip-tutorials). This lecture and other tutorial notebooks are indexed at the [QuTiP Tutorial webpage](https://qutip.org/tutorials.html).
 
@@ -410,7 +410,7 @@ psi0 = basis(2, 0)
 # list of times for which the solver should store the state vector
 tlist = np.linspace(0, 10, 100)
 
-result = mesolve(H, psi0, tlist, [], [])
+result = mesolve(H, psi0, tlist, [])
 ```
 
 ```python
@@ -458,12 +458,12 @@ axes.set_xlabel(r"$t$", fontsize=20)
 axes.set_ylabel(r"$\left<\sigma_z\right>$", fontsize=20);
 ```
 
-If we are only interested in expectation values, we could pass a list of operators to the `mesolve` function that we want expectation values for, and have the solver compute then and store the results in the `Odedata` class instance that it returns.
+If we are only interested in expectation values, we could pass a list of operators to the `mesolve` function that we want expectation values for, and have the solver compute then and store the results in the `Odedata` class instance that it returns.\
 
 For example, to request that the solver calculates the expectation values for the operators $\sigma_x, \sigma_y, \sigma_z$:
 
 ```python
-result = mesolve(H, psi0, tlist, [], [sigmax(), sigmay(), sigmaz()])
+result = mesolve(H, psi0, tlist, [], e_ops=[sigmax(), sigmay(), sigmaz()])
 ```
 
 Now the expectation values are available in `result.expect[0]`, `result.expect[1]`, and `result.expect[2]`:
diff --git a/tutorials-v5/lectures/Lecture-1-Jaynes-Cumming-model.md b/tutorials-v5/lectures/Lecture-1-Jaynes-Cumming-model.md
index 53b4c2c4..a097c6e4 100644
--- a/tutorials-v5/lectures/Lecture-1-Jaynes-Cumming-model.md
+++ b/tutorials-v5/lectures/Lecture-1-Jaynes-Cumming-model.md
@@ -17,7 +17,7 @@ jupyter:
 
 Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
 
-This lecture series was developed by J.R. Johannson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
+This lecture series was developed by J.R. Johansson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
 
 This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the [qutip-tutorials repository](https://github.com/qutip/qutip-tutorials). This lecture and other tutorial notebooks are indexed at the [QuTiP Tutorial webpage](https://qutip.org/tutorials.html).
 
@@ -69,7 +69,7 @@ tlist = np.linspace(0, 25, 101)
 ### Setup the operators, the Hamiltonian and initial state
 
 ```python
-# intial state
+# initial state
 psi0 = tensor(basis(N, 0), basis(2, 1))  # start with an excited atom
 
 # operators
@@ -111,7 +111,7 @@ if rate > 0.0:
 Here we evolve the system with the Lindblad master equation solver, and we request that the expectation values of the operators $a^\dagger a$ and $\sigma_+\sigma_-$ are returned by the solver by passing the list `[a.dag()*a, sm.dag()*sm]` as the fifth argument to the solver.
 
 ```python
-output = mesolve(H, psi0, tlist, c_ops, [a.dag() * a, sm.dag() * sm])
+output = mesolve(H, psi0, tlist, c_ops, e_ops=[a.dag() * a, sm.dag() * sm])
 ```
 
 ## Visualize the results
@@ -136,10 +136,10 @@ axes.set_title("Vacuum Rabi oscillations");
 
 In addition to the cavity's and atom's excitation probabilities, we may also be interested in for example the wigner function as a function of time. The Wigner function can give some valuable insight in the nature of the state of the resonators.
 
-To calculate the Wigner function in QuTiP, we first recalculte the evolution without specifying any expectation value operators, which will result in that the solver return a list of density matrices for the system for the given time coordinates.
+To calculate the Wigner function in QuTiP, we first recalculte the evolution without specifying any expectation value operators, which will result in that the solver returns a list of density matrices for the system for the given time coordinates.
 
 ```python
-output = mesolve(H, psi0, tlist, c_ops, [])
+output = mesolve(H, psi0, tlist, c_ops)
 ```
 
 Now, `output.states` contains a list of density matrices for the system for the time points specified in the list `tlist`:
@@ -167,7 +167,7 @@ For each of these points in time we need to:
 
  1. Find the system density matrix for the points in time that we are interested in.
  2. Trace out the atom and obtain the reduced density matrix for the cavity.
- 3. Calculate and visualize the Wigner function fo the reduced cavity density matrix.
+ 3. Calculate and visualize the Wigner function for the reduced cavity density matrix.
 
 ```python
 # find the indices of the density matrices for the times we are interested in
@@ -210,7 +210,7 @@ for idx, rho in enumerate(rho_list):
     axes[idx].set_title(r"$t = %.1f$" % tlist[t_idx][idx], fontsize=16)
 ```
 
-At $t =0$, the cavity is in it's ground state. At $t = 5, 15, 25$ it reaches it's maxium occupation in this Rabi-vacuum oscillation process. We can note that for $t=5$ and $t=15$ the Wigner function has negative values, indicating a truely quantum mechanical state. At $t=25$, however, the wigner function no longer has negative values and can therefore be considered a classical state.
+At $t =0$, the cavity is in its ground state. At $t = 5, 15, 25$ it reaches its maximum occupation in this Rabi-vacuum oscillation process. We can note that for $t=5$ and $t=15$ the Wigner function has negative values, indicating a truly quantum mechanical state. At $t=25$, however, the Wigner function no longer has negative values and can therefore be considered a classical state.
 
 Also, `qutip.anim_wigner` is useful to visualize the oscillation. It calculates the wigner function inside it so all you need to do is pass the reduced density matrices. You can pass options to it. For example, `projection='3d'` produces a 3d plot. You can see the oscillation easily with the 3d plot.
 
diff --git a/tutorials-v5/lectures/Lecture-10-cQED-dispersive-regime.md b/tutorials-v5/lectures/Lecture-10-cQED-dispersive-regime.md
index 78c1600e..d401b695 100644
--- a/tutorials-v5/lectures/Lecture-10-cQED-dispersive-regime.md
+++ b/tutorials-v5/lectures/Lecture-10-cQED-dispersive-regime.md
@@ -16,7 +16,7 @@ jupyter:
 
 Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
 
-This lecture series was developed by J.R. Johannson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
+This lecture series was developed by J.R. Johansson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
 
 This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the [qutip-tutorials repository](https://github.com/qutip/qutip-tutorials). This lecture and other tutorial notebooks are indexed at the [QuTiP Tutorial webpage](https://qutip.org/tutorials.html).
 
@@ -24,7 +24,7 @@ This is a slightly modified version of the lectures, to work with the current re
 import matplotlib as mpl
 import matplotlib.pyplot as plt
 import numpy as np
-from qutip import (SolverOptions, about, basis, coherent, correlation_2op_1t,
+from qutip import (about, basis, coherent, correlation_2op_1t,
                    destroy, expect, mesolve, ptrace, qeye, sigmax, sigmaz,
                    spectrum_correlation_fft, tensor, wigner)
 
@@ -115,7 +115,7 @@ tlist = np.linspace(0, 250, 1000)
 ```
 
 ```python
-res = mesolve(H, psi0, tlist, [], [], options=SolverOptions(nsteps=5000))
+res = mesolve(H, psi0, tlist, [], options={'nsteps': 5000})
 ```
 
 ### Excitation numbers
diff --git a/tutorials-v5/lectures/Lecture-11-Charge-Qubits.md b/tutorials-v5/lectures/Lecture-11-Charge-Qubits.md
index 8d3017cd..0440d87e 100644
--- a/tutorials-v5/lectures/Lecture-11-Charge-Qubits.md
+++ b/tutorials-v5/lectures/Lecture-11-Charge-Qubits.md
@@ -16,7 +16,7 @@ jupyter:
 
 Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
 
-This lecture series was developed by J.R. Johannson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
+This lecture series was developed by J.R. Johansson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
 
 This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the [qutip-tutorials repository](https://github.com/qutip/qutip-tutorials). This lecture and other tutorial notebooks are indexed at the [QuTiP Tutorial webpage](https://qutip.org/tutorials.html).
 
@@ -39,7 +39,7 @@ where $E_C$ is the charge energy, $E_J$ is the Josephson energy, and $\left| n\r
 
 #### References
 
- * [J. Koch et al, Phys. Rec. A 76, 042319 (2007)](https://journals.aps.org/pra/abstract/10.1103/PhysRevA.76.042319)
+ * [J. Koch et al, Phys. Rev. A 76, 042319 (2007)](https://journals.aps.org/pra/abstract/10.1103/PhysRevA.76.042319)
  * [Y.A. Pashkin et al, Quantum Inf Process 8, 55 (2009)](http://dx.doi.org/10.1007/s11128-009-0101-5)
 
 
@@ -176,7 +176,7 @@ energies = np.array([hamiltonian(Ec, Ej, N, ng).eigenenergies()
 plot_energies(ng_vec, energies, ymax=(50, 3));
 ```
 
-Note that the energy-level splitting is essentially independent of the gate bias $n_g$, at least for the lowest few states. This device insensitive to charge noise. But at the same time the two lowest energy states are no longer well separated from higher states (it has become more like an harmonic oscillator). But some anharmonicity still remains, and it can still be used as a qubit if the leakage of occupation probability of the higher states can be kept under control.
+Note that the energy-level splitting is essentially independent of the gate bias $n_g$, at least for the lowest few states. This device is insensitive to charge noise. But at the same time the two lowest energy states are no longer well separated from higher states (it has become more like a harmonic oscillator). But some anharmonicity still remains, and it can still be used as a qubit if the leakage of occupation probability to the higher states can be kept under control.
 
 
 ## Focus on the two lowest energy states
@@ -304,7 +304,7 @@ psi0 = psi_g
 
 ```python
 tlist = np.linspace(0.0, 100.0, 500)
-result = mesolve(Heff, psi0, tlist, [], [ket2dm(psi_e)], args=args)
+result = mesolve(Heff, psi0, tlist, [], e_ops=[ket2dm(psi_e)], args=args)
 ```
 
 ```python
@@ -362,7 +362,7 @@ psi_e = Qobj(psi_e.full()[keep_states, :])
 
 ```python
 tlist = np.linspace(0.0, 100.0, 500)
-result = mesolve(Heff, psi0, tlist, [], [ket2dm(psi_e)], args=args)
+result = mesolve(Heff, psi0, tlist, [], e_ops=[ket2dm(psi_e)], args=args)
 ```
 
 ```python
diff --git a/tutorials-v5/lectures/Lecture-12-Decay-into-a-squeezed-vacuum-field.md b/tutorials-v5/lectures/Lecture-12-Decay-into-a-squeezed-vacuum-field.md
index 5e669a09..0f89cb65 100644
--- a/tutorials-v5/lectures/Lecture-12-Decay-into-a-squeezed-vacuum-field.md
+++ b/tutorials-v5/lectures/Lecture-12-Decay-into-a-squeezed-vacuum-field.md
@@ -17,7 +17,7 @@ jupyter:
 
 Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
 
-This lecture series was developed by J.R. Johannson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
+This lecture series was developed by J.R. Johansson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
 
 This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the [qutip-tutorials repository](https://github.com/qutip/qutip-tutorials). This lecture and other tutorial notebooks are indexed at the [QuTiP Tutorial webpage](https://qutip.org/tutorials.html).
 
@@ -70,7 +70,7 @@ w_th = 0.0 * 2 * np.pi
 
 ```python
 # the number of average excitations in the
-# environment mode w0 at temperture w_th
+# environment mode w0 at temperature w_th
 Nth = n_thermal(w0, w_th)
 
 Nth
@@ -161,7 +161,7 @@ e_ops = [sigmax(), sigmay(), sigmaz()]
 ```
 
 ```python
-result1 = mesolve(L, psi0, tlist, [], e_ops)
+result1 = mesolve(L, psi0, tlist, [], e_ops=e_ops)
 ```
 
 ```python
@@ -205,7 +205,7 @@ c_ops = [np.sqrt(gamma0) *
 ```
 
 ```python
-result2 = mesolve(H, psi0, tlist, c_ops, e_ops)
+result2 = mesolve(H, psi0, tlist, c_ops, e_ops=e_ops)
 ```
 
 And we can verify that it indeed gives the same results:
@@ -292,7 +292,7 @@ c_ops = [np.sqrt(gamma0) *
 ```
 
 ```python
-result1 = mesolve(H, psi0, tlist, c_ops, e_ops)
+result1 = mesolve(H, psi0, tlist, c_ops, e_ops=e_ops)
 ```
 
 ```python
@@ -304,7 +304,7 @@ c_ops = [np.sqrt(gamma0) *
 ```
 
 ```python
-result2 = mesolve(H, psi0, tlist, c_ops, e_ops)
+result2 = mesolve(H, psi0, tlist, c_ops, e_ops=e_ops)
 ```
 
 ```python
diff --git a/tutorials-v5/lectures/Lecture-13-Resonance-flourescence.md b/tutorials-v5/lectures/Lecture-13-Resonance-flourescence.md
index fc269284..747474f0 100644
--- a/tutorials-v5/lectures/Lecture-13-Resonance-flourescence.md
+++ b/tutorials-v5/lectures/Lecture-13-Resonance-flourescence.md
@@ -12,12 +12,12 @@ jupyter:
     name: python3
 ---
 
-# Lecture 13 - Resonance flourescence
+# Lecture 13 - Resonance fluorescence
 
 
 Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
 
-This lecture series was developed by J.R. Johannson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
+This lecture series was developed by J.R. Johansson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
 
 This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the [qutip-tutorials repository](https://github.com/qutip/qutip-tutorials). This lecture and other tutorial notebooks are indexed at the [QuTiP Tutorial webpage](https://qutip.org/tutorials.html).
 
@@ -34,6 +34,10 @@ from qutip import (about, basis, correlation_2op_1t, mesolve, n_thermal, num,
 <!-- #region -->
 ## Introduction
 
+Resonance fluorescence, is the interaction between a two-level quantum system and a coherent driving field.
+
+The Hamiltonian for a two-level system under resonant driving, is given by:
+
 $\displaystyle H_L = -\frac{\Omega}{2}(\sigma_+ + \sigma_-)$
 
 
@@ -44,12 +48,18 @@ $\displaystyle \frac{d}{dt}\rho = -i[H_L, \rho] + \gamma_0(N+1)\left(\sigma_-\rh
 ### Problem definition in QuTiP
 
 ```python
+# Rabi frequency (strength of the driving field)
 Omega = 1.0 * 2 * np.pi
 ```
 
 ```python
+# decay rate of the two-level system
 gamma0 = 0.05
+
+# thermal frequency of the bath
 w_th = 0.0
+
+# thermal occupation number for the bath
 N = n_thermal(Omega, w_th)
 ```
 
@@ -75,7 +85,7 @@ psi0 = basis(2, 0)
 
 ```python
 tlist = np.linspace(0, 20 / (2 * np.pi), 200)
-result = mesolve(HL, psi0, tlist, c_ops, e_ops)
+result = mesolve(HL, psi0, tlist, c_ops, e_ops=e_ops)
 ```
 
 ```python
@@ -90,45 +100,56 @@ axes[0].plot(result.times, result.expect[2], "b",
 axes[0].legend()
 axes[0].set_ylim(-1, 1)
 
+axes[0].set_ylabel(r"$\langle\sigma_{x,y,z}\rangle$", fontsize=16)
+axes[0].set_title("Bloch Vector Components vs Time")
 
 axes[1].plot(result.times, result.expect[5], "b", label=r"$P_e$")
 
-# axes[1].set_ylabel(r'$\langle\sigma_z\rangle$', fontsize=16)
+axes[1].set_ylabel(r"$P_g$", fontsize=16)
 axes[1].set_xlabel("time", fontsize=16)
 axes[1].legend()
-axes[1].set_ylim(0, 1);
+axes[1].set_ylim(0, 1)
+axes[1].set_title("Ground State Population vs Time");
 ```
 
 ```python
 fig, ax = plt.subplots(1, 1, figsize=(12, 6), sharex=True)
 
-
 for idx, gamma0 in enumerate([0.1 * Omega, 0.5 * Omega, 1.0 * Omega]):
 
     HL, c_ops = system_spec(Omega, gamma0, N)
-    result = mesolve(HL, psi0, tlist, c_ops, e_ops)
+    result = mesolve(HL, psi0, tlist, c_ops, e_ops=e_ops)
 
-    ax.plot(result.times, result.expect[5], "b",
-            label=r"$\langle\sigma_z\rangle$")
+    ax.plot(
+        result.times, result.expect[5], "b",
+        label=fr"$P_g$ ($\gamma_0={gamma0:.2f}$)"
+    )
 
-ax.set_ylim(0, 1);
+ax.set_ylim(0, 1)
+ax.set_xlabel("time", fontsize=16)
+ax.set_ylabel(r"$P_g$", fontsize=16)
+ax.set_title("Ground State Population for Different Decay Rates")
+ax.legend();
 ```
 
 ```python
 fig, ax = plt.subplots(1, 1, figsize=(12, 6), sharex=True)
 
-
 for idx, gamma0 in enumerate([0.1 * Omega, 0.5 * Omega, 1.0 * Omega]):
 
     HL, c_ops = system_spec(Omega, gamma0, N)
-    result = mesolve(HL, psi0, tlist, c_ops, e_ops)
+    result = mesolve(HL, psi0, tlist, c_ops, e_ops=e_ops)
 
     ax.plot(
         result.times, np.imag(result.expect[4]),
-        label=r"im $\langle\sigma_+\rangle$"
+        label=fr"Im $\langle\sigma_+\rangle$ ($\gamma_0={gamma0:.2f}$)"
     )
 
-ax.set_ylim(-0.5, 0.5);
+ax.set_ylim(-0.5, 0.5)
+ax.set_xlabel("time", fontsize=16)
+ax.set_ylabel(r"Im $\langle\sigma_+\rangle$", fontsize=16)
+ax.set_title(r"Imaginary Part of $\langle\sigma_+\rangle$ for Different Decay Rates")
+ax.legend();
 ```
 
 ```python
@@ -142,12 +163,21 @@ for idx, gamma0 in enumerate([2 * Omega, 0.5 * Omega, 0.25 * Omega]):
     corr_vec = correlation_2op_1t(HL, None, taulist, c_ops, sigmap(), sigmam())
     w, S = spectrum_correlation_fft(taulist, corr_vec)
 
-    axes[0].plot(taulist, corr_vec, label=r"$<\sigma_+(\tau)\sigma_-(0)>$")
-    axes[1].plot(-w / (gamma0), S, "b", label=r"$S(\omega)$")
-    axes[1].plot(w / (gamma0), S, "b", label=r"$S(\omega)$")
+    axes[0].plot(taulist, corr_vec, label=fr"$\gamma_0={gamma0:.2f}$")
+    axes[1].plot(-w / (gamma0), S, "b", label=fr"$\gamma_0={gamma0:.2f}$")
+    axes[1].plot(w / (gamma0), S, "b")
 
 axes[0].set_xlim(0, 10)
-axes[1].set_xlim(-5, 5);
+axes[0].set_xlabel(r"$\tau$", fontsize=16)
+axes[0].set_ylabel(r"$\langle\sigma_+(\tau)\sigma_-(0)\rangle$", fontsize=16)
+axes[0].set_title("Two-Time Correlation Function")
+axes[0].legend()
+
+axes[1].set_xlim(-5, 5)
+axes[1].set_xlabel(r"$\omega/\gamma_0$", fontsize=16)
+axes[1].set_ylabel(r"$S(\omega)$", fontsize=16)
+axes[1].set_title("Resonance Fluorescence Spectrum")
+axes[1].legend();
 ```
 
 ### Software versions
diff --git a/tutorials-v5/lectures/Lecture-14-Kerr-nonlinearities.md b/tutorials-v5/lectures/Lecture-14-Kerr-nonlinearities.md
index b6bb4288..d03b0ed7 100644
--- a/tutorials-v5/lectures/Lecture-14-Kerr-nonlinearities.md
+++ b/tutorials-v5/lectures/Lecture-14-Kerr-nonlinearities.md
@@ -16,7 +16,7 @@ jupyter:
 
 Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
 
-This lecture series was developed by J.R. Johannson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
+This lecture series was developed by J.R. Johansson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
 
 This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the [qutip-tutorials repository](https://github.com/qutip/qutip-tutorials). This lecture and other tutorial notebooks are indexed at the [QuTiP Tutorial webpage](https://qutip.org/tutorials.html).
 
@@ -40,7 +40,7 @@ The Kerr effect describes a self-interaction electromagnetic quantum field which
 
 $\displaystyle H = \frac{1}{2}\chi (a^\dagger)^2a^2$
 
-where $\chi$ is related to the third-order nonlinear suseptibility. The Kerr effect is one of the typical nonlinearities that can occur in quantum optics due to a nonlinear medium.
+where $\chi$ is related to the third-order nonlinear susceptibility. The Kerr effect is one of the typical nonlinearities that can occur in quantum optics due to a nonlinear medium.
 
 In this notebook we'll see how to setup the model in QuTiP and look at some interesting properties of the states that evolve according to this Hamiltonian.
 
@@ -123,7 +123,7 @@ def plot_wigner(rho, fig=None, ax=None):
         W,
         100,
         norm=mpl.colors.Normalize(-wlim, wlim),
-        cmap=mpl.cm.get_cmap("RdBu"),
+        cmap=mpl.colormaps["RdBu"],
     )
     ax.set_xlabel(r"$x_1$", fontsize=16)
     ax.set_ylabel(r"$x_2$", fontsize=16)
@@ -148,7 +148,7 @@ def plot_fock_distribution_vs_time(tlist, states, fig=None, ax=None):
         Y,
         Z.T,
         norm=mpl.colors.Normalize(0, 0.5),
-        cmap=mpl.cm.get_cmap("Reds"),
+        cmap=mpl.colormaps["RdBu"],
         edgecolors="k",
     )
     ax.set_xlabel(r"$N$", fontsize=16)
@@ -184,7 +184,7 @@ psi0 = coherent(N, 2.0)
 # and evolve the state under the influence of the hamiltonian.
 # by passing an empty list as expecation value operators argument,
 # we get the full state of the system in result.states
-result = mesolve(H, psi0, tlist, [], [])
+result = mesolve(H, psi0, tlist, [])
 ```
 
 First, let's look at how the expecation values and variances of the photon number operator $n$ and the $x$ and $p$ quadratures evolve in time:
diff --git a/tutorials-v5/lectures/Lecture-15-Nonclassically-driven-atoms.md b/tutorials-v5/lectures/Lecture-15-Nonclassically-driven-atoms.md
index 8ca0956d..0370923f 100644
--- a/tutorials-v5/lectures/Lecture-15-Nonclassically-driven-atoms.md
+++ b/tutorials-v5/lectures/Lecture-15-Nonclassically-driven-atoms.md
@@ -17,7 +17,7 @@ jupyter:
 
 Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
 
-This lecture series was developed by J.R. Johannson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
+This lecture series was developed by J.R. Johansson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
 
 This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the [qutip-tutorials repository](https://github.com/qutip/qutip-tutorials). This lecture and other tutorial notebooks are indexed at the [QuTiP Tutorial webpage](https://qutip.org/tutorials.html).
 
diff --git a/tutorials-v5/lectures/Lecture-16-Gallery-of-Wigner-functions.md b/tutorials-v5/lectures/Lecture-16-Gallery-of-Wigner-functions.md
index 997adce3..ff186a89 100644
--- a/tutorials-v5/lectures/Lecture-16-Gallery-of-Wigner-functions.md
+++ b/tutorials-v5/lectures/Lecture-16-Gallery-of-Wigner-functions.md
@@ -17,7 +17,7 @@ jupyter:
 
 Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
 
-This lecture series was developed by J.R. Johannson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
+This lecture series was developed by J.R. Johansson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
 
 This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the [qutip-tutorials repository](https://github.com/qutip/qutip-tutorials). This lecture and other tutorial notebooks are indexed at the [QuTiP Tutorial webpage](https://qutip.org/tutorials.html).
 
@@ -132,10 +132,6 @@ plot_wigner_2d_3d(psi)
 
 ## Fock states: $\left|n\right>$
 
-```python
-
-```
-
 ```python
 for n in range(6):
     psi = basis(N, n)
diff --git a/tutorials-v5/lectures/Lecture-2A-Cavity-Qubit-Gates.md b/tutorials-v5/lectures/Lecture-2A-Cavity-Qubit-Gates.md
index 0b60a474..164f7fae 100644
--- a/tutorials-v5/lectures/Lecture-2A-Cavity-Qubit-Gates.md
+++ b/tutorials-v5/lectures/Lecture-2A-Cavity-Qubit-Gates.md
@@ -16,7 +16,7 @@ jupyter:
 
 Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
 
-This lecture series was developed by J.R. Johannson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
+This lecture series was developed by J.R. Johansson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
 
 This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the [qutip-tutorials repository](https://github.com/qutip/qutip-tutorials). This lecture and other tutorial notebooks are indexed at the [QuTiP Tutorial webpage](https://qutip.org/tutorials.html).
 
@@ -69,7 +69,7 @@ sm1 = tensor(qeye(N), destroy(2), qeye(2))
 sz1 = tensor(qeye(N), sigmaz(), qeye(2))
 n1 = sm1.dag() * sm1
 
-# oeprators for qubit 2
+# operators for qubit 2
 sm2 = tensor(qeye(N), qeye(2), destroy(2))
 sz2 = tensor(qeye(N), qeye(2), sigmaz())
 n2 = sm2.dag() * sm2
@@ -139,7 +139,7 @@ H_t = [[Hc, wc_t], [H1, w1_t], [H2, w2_t], Hc1 + Hc2]
 ### Evolve the system
 
 ```python
-res = mesolve(H_t, psi0, tlist, [], [])
+res = mesolve(H_t, psi0, tlist, [])
 ```
 
 ### Plot the results
@@ -239,7 +239,7 @@ c_ops = [np.sqrt(kappa) * a, np.sqrt(gamma1) * sm1, np.sqrt(gamma2) * sm2]
 ### Evolve the system
 
 ```python
-res = mesolve(H_t, psi0, tlist, c_ops, [])
+res = mesolve(H_t, psi0, tlist, c_ops)
 ```
 
 ### Plot the results
@@ -323,7 +323,7 @@ fig.tight_layout()
 ### Evolve the system
 
 ```python
-res = mesolve(H_t, psi0, tlist, [], [])
+res = mesolve(H_t, psi0, tlist, [])
 ```
 
 ### Plot the results
@@ -408,7 +408,7 @@ fig.tight_layout()
 ### Evolve the system
 
 ```python
-res = mesolve(H_t, psi0, tlist, [], [])
+res = mesolve(H_t, psi0, tlist, [])
 ```
 
 ### Plot the results
@@ -488,7 +488,7 @@ c_ops = [np.sqrt(kappa) * a, np.sqrt(gamma1) * sm1, np.sqrt(gamma2) * sm2]
 ### Evolve the system
 
 ```python
-res = mesolve(H_t, psi0, tlist, c_ops, [])
+res = mesolve(H_t, psi0, tlist, c_ops)
 ```
 
 ### Plot results
@@ -574,7 +574,7 @@ H_t = [[Hc, wc_t], H1 * w1 + H2 * w2 + Hc1 + Hc2]
 ### Evolve the system
 
 ```python
-res = mesolve(H_t, psi0, tlist, c_ops, [])
+res = mesolve(H_t, psi0, tlist, c_ops)
 ```
 
 ### Plot the results
diff --git a/tutorials-v5/lectures/Lecture-2B-Single-Atom-Lasing.md b/tutorials-v5/lectures/Lecture-2B-Single-Atom-Lasing.md
index 56191917..9b86899c 100644
--- a/tutorials-v5/lectures/Lecture-2B-Single-Atom-Lasing.md
+++ b/tutorials-v5/lectures/Lecture-2B-Single-Atom-Lasing.md
@@ -17,7 +17,7 @@ jupyter:
 
 Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
 
-This lecture series was developed by J.R. Johannson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
+This lecture series was developed by J.R. Johansson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
 
 This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the [qutip-tutorials repository](https://github.com/qutip/qutip-tutorials). This lecture and other tutorial notebooks are indexed at the [QuTiP Tutorial webpage](https://qutip.org/tutorials.html).
 
@@ -28,7 +28,7 @@ import matplotlib as mpl
 import matplotlib.pyplot as plt
 import numpy as np
 from IPython.display import Image
-from qutip import (SolverOptions, about, basis, destroy, expect, mesolve,
+from qutip import (about, basis, destroy, expect, mesolve,
                    ptrace, qeye, sigmax, steadystate, tensor, wigner)
 
 %matplotlib inline
@@ -90,7 +90,7 @@ tlist = np.linspace(0, 150, 101)
 ### Setup the operators, the Hamiltonian and initial state
 
 ```python
-# intial state
+# initial state
 psi0 = tensor(basis(N, 0), basis(2, 0))  # start without excitations
 
 # operators
@@ -130,8 +130,8 @@ if rate > 0.0:
 Here we evolve the system with the Lindblad master equation solver, and we request that the expectation values of the operators $a^\dagger a$ and $\sigma_+\sigma_-$ are returned by the solver by passing the list `[a.dag()*a, sm.dag()*sm]` as the fifth argument to the solver.
 
 ```python
-opt = SolverOptions(nsteps=2000)  # allow extra time-steps
-output = mesolve(H, psi0, tlist, c_ops, [a.dag() * a, sm.dag() * sm],
+opt = {'nsteps': 2000}  # allow extra time-steps
+output = mesolve(H, psi0, tlist, c_ops, e_ops=[a.dag() * a, sm.dag() * sm],
                  options=opt)
 ```
 
@@ -191,8 +191,8 @@ axes[0].set_ylabel("Occupation probability", fontsize=18);
 
 ```python
 tlist = np.linspace(0, 25, 5)
-output = mesolve(H, psi0, tlist, c_ops, [],
-                 options=SolverOptions(nsteps=5000))
+output = mesolve(H, psi0, tlist, c_ops,
+                 options={'nsteps': 5000})
 ```
 
 ```python
diff --git a/tutorials-v5/lectures/Lecture-3A-Dicke-model.md b/tutorials-v5/lectures/Lecture-3A-Dicke-model.md
index 3aaae2ab..1049fce4 100644
--- a/tutorials-v5/lectures/Lecture-3A-Dicke-model.md
+++ b/tutorials-v5/lectures/Lecture-3A-Dicke-model.md
@@ -16,7 +16,7 @@ jupyter:
 
 Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
 
-This lecture series was developed by J.R. Johannson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
+This lecture series was developed by J.R. Johansson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
 
 This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the [qutip-tutorials repository](https://github.com/qutip/qutip-tutorials). This lecture and other tutorial notebooks are indexed at the [QuTiP Tutorial webpage](https://qutip.org/tutorials.html).
 
@@ -201,7 +201,7 @@ fig.tight_layout()
 * [Lambert et al., Phys. Rev. Lett. 92, 073602 (2004)](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.92.073602).
 
 ```python
-def calulcate_entropy(M, N, g_vec):
+def calculate_entropy(M, N, g_vec):
 
     j = N / 2.0
     n = 2 * j + 1
@@ -241,7 +241,7 @@ fig, axes = plt.subplots(1, 1, figsize=(12, 6))
 
 for NN in N_vec:
 
-    entropy_cavity, entropy_spin = calulcate_entropy(MM, NN, g_vec)
+    entropy_cavity, entropy_spin = calculate_entropy(MM, NN, g_vec)
 
     axes.plot(g_vec, entropy_cavity, "b", label="N = %d" % NN)
     axes.plot(g_vec, entropy_spin, "r--")
diff --git a/tutorials-v5/lectures/Lecture-3B-Jaynes-Cumming-model-with-ultrastrong-coupling.md b/tutorials-v5/lectures/Lecture-3B-Jaynes-Cumming-model-with-ultrastrong-coupling.md
index cce8fc10..1d31817a 100644
--- a/tutorials-v5/lectures/Lecture-3B-Jaynes-Cumming-model-with-ultrastrong-coupling.md
+++ b/tutorials-v5/lectures/Lecture-3B-Jaynes-Cumming-model-with-ultrastrong-coupling.md
@@ -118,7 +118,7 @@ na_expt = expect(na, psi_list)  # qubit  occupation probability
 nc_expt = expect(nc, psi_list)  # cavity occupation probability
 ```
 
-Plot the ground state occupation probabilities of the cavity and the atom as a function of coupling strenght. Note that for large coupling strength (the ultrastrong coupling regime, where $g > \omega_a,\omega_c$), the ground state has both photonic and atomic excitations.
+Plot the ground state occupation probabilities of the cavity and the atom as a function of coupling strength. Note that for large coupling strength (the ultrastrong coupling regime, where $g > \omega_a,\omega_c$), the ground state has both photonic and atomic excitations.
 
 ```python
 fig, axes = plt.subplots(1, 1, sharex=True, figsize=(8, 4))
@@ -126,7 +126,7 @@ fig, axes = plt.subplots(1, 1, sharex=True, figsize=(8, 4))
 axes.plot(g_vec / (2 * np.pi), nc_expt, "r", linewidth=2, label="cavity")
 axes.plot(g_vec / (2 * np.pi), na_expt, "b", linewidth=2, label="atom")
 axes.set_ylabel("Occupation probability", fontsize=16)
-axes.set_xlabel("coupling strenght", fontsize=16)
+axes.set_xlabel("coupling strength", fontsize=16)
 axes.legend(loc=0)
 
 fig.tight_layout()
@@ -206,7 +206,7 @@ psi0 = tensor(basis(N, 1), basis(2, 0))
 
 ```python
 tlist = np.linspace(0, 20, 1000)
-output = mesolve(H, psi0, tlist, [], [a.dag() * a, sm.dag() * sm])
+output = mesolve(H, psi0, tlist, [], e_ops=[a.dag() * a, sm.dag() * sm])
 ```
 
 ```python
@@ -223,7 +223,7 @@ fig.tight_layout()
 
 ```python
 tlist = np.linspace(0, 0.35, 8)
-output = mesolve(H, psi0, tlist, [], [])
+output = mesolve(H, psi0, tlist, [])
 ```
 
 ```python
@@ -269,7 +269,7 @@ kappa = 0.25
 ```python
 tlist = np.linspace(0, 20, 1000)
 output = mesolve(H, psi0, tlist, [np.sqrt(kappa) * a],
-                 [a.dag() * a, sm.dag() * sm])
+                 e_ops=[a.dag() * a, sm.dag() * sm])
 ```
 
 ```python
@@ -281,7 +281,7 @@ axes.legend(loc=0);
 
 ```python
 tlist = np.linspace(0, 10, 8)
-output = mesolve(H, psi0, tlist, [np.sqrt(kappa) * a], [])
+output = mesolve(H, psi0, tlist, [np.sqrt(kappa) * a])
 ```
 
 ```python
@@ -323,7 +323,7 @@ tlist = np.linspace(0, 30, 50)
 
 psi0 = H.groundstate()[1]
 
-output = mesolve(H, psi0, tlist, [np.sqrt(kappa) * a], [])
+output = mesolve(H, psi0, tlist, [np.sqrt(kappa) * a])
 ```
 
 ```python
diff --git a/tutorials-v5/lectures/Lecture-4-Correlation-Functions.md b/tutorials-v5/lectures/Lecture-4-Correlation-Functions.md
index 3d4990bc..1cd86510 100644
--- a/tutorials-v5/lectures/Lecture-4-Correlation-Functions.md
+++ b/tutorials-v5/lectures/Lecture-4-Correlation-Functions.md
@@ -16,7 +16,7 @@ jupyter:
 
 Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
 
-This lecture series was developed by J.R. Johannson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
+This lecture series was developed by J.R. Johansson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
 
 This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the [qutip-tutorials repository](https://github.com/qutip/qutip-tutorials). This lecture and other tutorial notebooks are indexed at the [QuTiP Tutorial webpage](https://qutip.org/tutorials.html).
 
@@ -61,7 +61,7 @@ c_ops = [np.sqrt(G1 * (1 + n_th)) * a, np.sqrt(G1 * n_th) * a.dag()]
 rho0 = coherent_dm(N, 2.0)
 
 # first calculate the occupation number as a function of time
-n = mesolve(H, rho0, taulist, c_ops, [a.dag() * a]).expect[0]
+n = mesolve(H, rho0, taulist, c_ops, e_ops=[a.dag() * a]).expect[0]
 n = np.array(n)
 
 # calculate the correlation function G1 and normalize with n to obtain g1
@@ -112,7 +112,7 @@ def correlation_ss_gtt(H, tlist, c_ops, a_op, b_op, c_op, d_op, rho0=None):
         rho0 = steadystate(H, c_ops)
 
     return mesolve(H, d_op * rho0 * a_op, tlist, c_ops,
-                   [b_op * c_op]).expect[0]
+                   e_ops=[b_op * c_op]).expect[0]
 ```
 
 ```python
diff --git a/tutorials-v5/lectures/Lecture-5-Parametric-Amplifier.md b/tutorials-v5/lectures/Lecture-5-Parametric-Amplifier.md
index e9af16e4..9ba86be0 100644
--- a/tutorials-v5/lectures/Lecture-5-Parametric-Amplifier.md
+++ b/tutorials-v5/lectures/Lecture-5-Parametric-Amplifier.md
@@ -16,7 +16,7 @@ jupyter:
 
 Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
 
-This lecture series was developed by J.R. Johannson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
+This lecture series was developed by J.R. Johansson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
 
 This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the [qutip-tutorials repository](https://github.com/qutip/qutip-tutorials). This lecture and other tutorial notebooks are indexed at the [QuTiP Tutorial webpage](https://qutip.org/tutorials.html).
 
@@ -76,7 +76,7 @@ e_ops = []
 ```
 
 ```python
-output = mesolve(H, psi0, tlist, c_ops, e_ops)
+output = mesolve(H, psi0, tlist, c_ops, e_ops=e_ops)
 output
 ```
 
@@ -95,7 +95,7 @@ for idx, psi in enumerate(output.states):
     nb_e[idx] = expect(nb, psi)
     nb_s[idx] = expect(nb * nb, psi)
 
-# substract the average squared to obtain variances
+# subtract the average squared to obtain variances
 na_s = na_s - na_e**2
 nb_s = nb_s - nb_e**2
 ```
@@ -190,10 +190,16 @@ for idx, psi in enumerate(output.states):
     cs_lhs[idx] = expect(ad_a_bd_b, psi)
     cs_rhs[idx] = expect(ad_ad_a_a, psi)
 
-# normalize the correlation functions
-g2_1 = g2_1 / (na_e**2)
-g2_2 = g2_2 / (nb_e**2)
-g2_12 = g2_12 / (na_e * nb_e)
+# normalize setting inf to nan
+def safe_divide(a, b):
+    with np.errstate(divide='ignore', invalid='ignore'):
+        result = np.true_divide(a, b)
+        result[~np.isfinite(result)] = np.nan
+    return result
+
+g2_1 = safe_divide(g2_1, na_e**2)
+g2_2 = safe_divide(g2_2, nb_e**2)
+g2_12 = safe_divide(g2_12, na_e * nb_e)
 ```
 
 ### Second-order coherence functions: Cauchy-Schwarz inequality
@@ -248,7 +254,9 @@ line1 = axes[0].plot(tlist, cs_lhs, "b", tlist, cs_rhs, "r", linewidth=2)
 axes[0].set_xlabel("$t$", fontsize=18)
 axes[0].set_title(r"Cauchy-Schwarz inequality", fontsize=18)
 
-line1 = axes[1].plot(tlist, cs_lhs / (cs_rhs), "k", linewidth=2)
+cs_ratio = safe_divide(cs_lhs, cs_rhs)
+
+line1 = axes[1].plot(tlist, cs_ratio, "k", linewidth=2)
 axes[1].set_xlabel("$t$", fontsize=18)
 axes[1].set_title(r"Cauchy-Schwarz ratio inequality", fontsize=18)
 
@@ -416,12 +424,11 @@ def plot_covariance_matrix(V, ax):
 
     ax.view_init(azim=-40, elev=60)
     ax.bar3d(xpos, ypos, zpos, dx, dy, dz, color=colors)
-    ax.axes.xaxis.set_major_locator(plt.IndexLocator(1, -0.5))
-    ax.axes.yaxis.set_major_locator(plt.IndexLocator(1, -0.5))
-    ax.axes.xaxis.set_ticklabels(("$q_-$", "$p_-$", "$q_+$", "$p_+$"),
-                                 fontsize=12)
-    ax.axes.yaxis.set_ticklabels(("$q_-$", "$p_-$", "$q_+$", "$p_+$"),
-                                 fontsize=12)
+    # Set tick locations before setting tick labels
+    ax.axes.xaxis.set_ticks([0, 1, 2, 3])
+    ax.axes.yaxis.set_ticks([0, 1, 2, 3])
+    ax.axes.xaxis.set_ticklabels(("$q_-$", "$p_-$", "$q_+$", "$p_+$"), fontsize=12)
+    ax.axes.yaxis.set_ticklabels(("$q_-$", "$p_-$", "$q_+$", "$p_+$"), fontsize=12)
 ```
 
 ```python
@@ -440,7 +447,7 @@ for idx, t_idx in enumerate(t_idx_vec):
 
     plot_covariance_matrix(V, axes[idx])
 
-fig.tight_layout()
+fig.subplots_adjust(left=0.15, right=0.85, top=0.9, bottom=0.1)
 ```
 
 ```python
diff --git a/tutorials-v5/lectures/Lecture-6-Quantum-Monte-Carlo-Trajectories.md b/tutorials-v5/lectures/Lecture-6-Quantum-Monte-Carlo-Trajectories.md
index fd5c3a57..4f35f622 100644
--- a/tutorials-v5/lectures/Lecture-6-Quantum-Monte-Carlo-Trajectories.md
+++ b/tutorials-v5/lectures/Lecture-6-Quantum-Monte-Carlo-Trajectories.md
@@ -17,7 +17,7 @@ jupyter:
 
 Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
 
-This lecture series was developed by J.R. Johannson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
+This lecture series was developed by J.R. Johansson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
 
 This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the [qutip-tutorials repository](https://github.com/qutip/qutip-tutorials). This lecture and other tutorial notebooks are indexed at the [QuTiP Tutorial webpage](https://qutip.org/tutorials.html).
 
@@ -128,13 +128,13 @@ The expectation values of $a^\dagger a$ are now available in array ``mc.expect[i
 
 ## Lindblad master-equation simulation and steady state
 
-For comparison with the averages of single quantum trajectories provided by the Monte-Carlo solver we here also calculate the dynamics of the Lindblad master equation, which should agree with the Monte-Carlo simultions for infinite number of trajectories.
+For comparison with the averages of single quantum trajectories provided by the Monte-Carlo solver we here also calculate the dynamics of the Lindblad master equation, which should agree with the Monte-Carlo simulations for infinite number of trajectories.
 
 ```python
 # run master equation to get ensemble average expectation values
 me = mesolve(H, psi0, tlist, c_op_list, e_ops=[a.dag() * a])
 
-# calulate final state using steadystate solver
+# calculate final state using steadystate solver
 final_state = steadystate(H, c_op_list)  # find steady-state
 # find expectation value for particle number
 fexpt = expect(a.dag() * a, final_state)
diff --git a/tutorials-v5/lectures/Lecture-7-iSWAP-gate.md b/tutorials-v5/lectures/Lecture-7-iSWAP-gate.md
index ece8e269..e873a654 100644
--- a/tutorials-v5/lectures/Lecture-7-iSWAP-gate.md
+++ b/tutorials-v5/lectures/Lecture-7-iSWAP-gate.md
@@ -16,7 +16,7 @@ jupyter:
 
 Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
 
-This lecture series was developed by J.R. Johannson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
+This lecture series was developed by J.R. Johansson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
 
 This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the [qutip-tutorials repository](https://github.com/qutip/qutip-tutorials). This lecture and other tutorial notebooks are indexed at the [QuTiP Tutorial webpage](https://qutip.org/tutorials.html).
 
@@ -37,7 +37,7 @@ $\displaystyle H = g \left(\sigma_x\otimes\sigma_x + \sigma_y\otimes\sigma_y\rig
 
 where $g$ is the coupling strength. Under ideal conditions this coupling realizes the $i$-SWAP gate between the two qubit states. 
 
-Here we will solve for the dynamics of the two qubits subject to this Hamiltonian, and look at the deterioating effects of adding decoherence. We will use process tomography to visualize the gate.
+Here we will solve for the dynamics of the two qubits subject to this Hamiltonian, and look at the deteriorating effects of adding decoherence. We will use process tomography to visualize the gate.
 
 
 ### Parameters
diff --git a/tutorials-v5/lectures/Lecture-8-Adiabatic-quantum-computing.md b/tutorials-v5/lectures/Lecture-8-Adiabatic-quantum-computing.md
index eba45166..ec57478b 100644
--- a/tutorials-v5/lectures/Lecture-8-Adiabatic-quantum-computing.md
+++ b/tutorials-v5/lectures/Lecture-8-Adiabatic-quantum-computing.md
@@ -16,20 +16,20 @@ jupyter:
 
 Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
 
-This lecture series was developed by J.R. Johannson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
+This lecture series was developed by J.R. Johansson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
 
 This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the [qutip-tutorials repository](https://github.com/qutip/qutip-tutorials). This lecture and other tutorial notebooks are indexed at the [QuTiP Tutorial webpage](https://qutip.org/tutorials.html).
 
 
 ## Introduction
 
-In adiabatic quantum computing, an easy to prepare ground state of a Hamiltonian $H_0$ is prepared, and then the Hamiltonian is gradually transformed into $H_1$, which is constructed in such a way that the groundstate of $H_1$ encodes the solution to a difficult problem. The transformation of $H_0$ to $H_1$ can for example be written on the form
+In adiabatic quantum computing, an easy to prepare ground state of a Hamiltonian $H_0$ is prepared, and then the Hamiltonian is gradually transformed into $H_1$, which is constructed in such a way that the ground state of $H_1$ encodes the solution to a difficult problem. The transformation of $H_0$ to $H_1$ can for example be written on the form
 
 $\displaystyle H(t) = \lambda(t) H_0 + (1 - \lambda(t)) H_1$
 
-where $\lambda(t)$ is a function that goes from goes from $0$ to $1$ when $t$ goes from $0$ to $t_{\rm final}$.
+where $\lambda(t)$ is a function that goes from $0$ to $1$ when $t$ goes from $0$ to $t_{\rm final}$.
 
-If this gradual tranformation is slow enough (satisfying the adiabicity critera), the evolution of the system will remain in its ground state.
+If this gradual transformation is slow enough (satisfying the adiabaticity criteria), the evolution of the system will remain in its ground state.
 
 If the Hamiltonian is transformed from $H_0$ to $H_1$ too quickly, the system will get excited from the ground state the adiabatic computing algorithm fails.
 
diff --git a/tutorials-v5/lectures/Lecture-9-Squeezed-states-of-harmonic-oscillator.md b/tutorials-v5/lectures/Lecture-9-Squeezed-states-of-harmonic-oscillator.md
index 41767993..aa74b4f7 100644
--- a/tutorials-v5/lectures/Lecture-9-Squeezed-states-of-harmonic-oscillator.md
+++ b/tutorials-v5/lectures/Lecture-9-Squeezed-states-of-harmonic-oscillator.md
@@ -16,7 +16,7 @@ jupyter:
 
 Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
 
-This lecture series was developed by J.R. Johannson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
+This lecture series was developed by J.R. Johansson. The original lecture notebooks are available [here](https://github.com/jrjohansson/qutip-lectures).
 
 This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the [qutip-tutorials repository](https://github.com/qutip/qutip-tutorials). This lecture and other tutorial notebooks are indexed at the [QuTiP Tutorial webpage](https://qutip.org/tutorials.html).
 
@@ -35,19 +35,19 @@ from qutip import (about, basis, coherent, destroy, displace, expect, mesolve,
 
 ## Introduction
 
-In quantum mechanics, each measurement of an observable (corresponding to an Hermitian operator) results in stochastic outcome that follows some probability distribution. The expectation value of the operator is the average of many measurement outcomes, and the standard deviation of the operator describes the uncertainty in the outcomes.
+In quantum mechanics, each measurement of an observable (corresponding to a Hermitian operator) results in stochastic outcome that follows some probability distribution. The expectation value of the operator is the average of many measurement outcomes, and the standard deviation of the operator describes the uncertainty in the outcomes.
 
-This uncertainty is intrinsic in quantum mechanics, and cannot be eliminated. The Heisenberg uncertainty principle describes the minumum uncertainly for pairs of noncommuting operators. For example, the operators such $x$ and $p$, which satisfy the commutation relation $[x, p] = i\hbar$, must always satisfy $(\Delta x) (\Delta p) >= \hbar/2$ .
+This uncertainty is intrinsic in quantum mechanics, and cannot be eliminated. The Heisenberg uncertainty principle describes the minimum uncertainty for pairs of noncommuting operators. For example, the operators such $x$ and $p$, which satisfy the commutation relation $[x, p] = i\hbar$, must always satisfy $(\Delta x) (\Delta p) >= \hbar/2$ .
 
 A state that satisfies
 
 $(\Delta x) (\Delta p) = \hbar/2$
 
-is called a minimum uncertainty state, and a state for which, for example,
+is called a minimum uncertainty state, and a state for which, for example, 
 
-$(\Delta x)^2 < \hbar/2$
+$(\Delta x)^2 < \hbar/2$ 
 
-is called a squeezed state. Note that in this case $(\Delta p)^2$ must be larger than $\hbar/2(\Delta x)^2$ for the Heisenberg relation to hold. Squeezing a quantum state so that the variance of one operator $x$ is reduced below the minimum uncertainty limit therefore necessarily amplify the variance of operators that do not commute with $x$, such as $p$.
+is called a squeezed state. Note that in this case $(\Delta p)^2$ must be larger than $\hbar/2(\Delta x)^2$ for the Heisenberg relation to hold. Squeezing a quantum state so that the variance of one operator $x$ is reduced below the minimum uncertainty limit therefore necessarily amplifies the variance of operators that do not commute with $x$, such as $p$.
 
 For harmonic modes, squeezing of $x$ or $p$ is called quadrature squeezing, and it is probably the most common form of squeezing.
 
@@ -130,7 +130,7 @@ psi0 = coherent(N, 2.0)
 ```
 
 ```python
-result = mesolve(H, psi0, tlist, c_ops, [])
+result = mesolve(H, psi0, tlist, c_ops)
 ```
 
 ```python
@@ -168,7 +168,7 @@ psi0 = squeeze(N, 1.0) * basis(N, 0)
 ```
 
 ```python
-result = mesolve(H, psi0, tlist, c_ops, [])
+result = mesolve(H, psi0, tlist, c_ops)
 ```
 
 ```python
@@ -208,7 +208,7 @@ psi0 = (
 ```
 
 ```python
-result = mesolve(H, psi0, tlist, c_ops, [])
+result = mesolve(H, psi0, tlist, c_ops)
 ```
 
 ```python