diff --git a/.github/workflows/Invalidations.yml b/.github/workflows/Invalidations.yml
deleted file mode 100644
index 4616f1cb..00000000
--- a/.github/workflows/Invalidations.yml
+++ /dev/null
@@ -1,40 +0,0 @@
-name: Invalidations
-
-on:
-  pull_request:
-
-concurrency:
-  # Skip intermediate builds: always.
-  # Cancel intermediate builds: always.
-  group: ${{ github.workflow }}-${{ github.ref }}
-  cancel-in-progress: true
-
-jobs:
-  evaluate:
-    # Only run on PRs to the default branch.
-    # In the PR trigger above branches can be specified only explicitly whereas this check should work for master, main, or any other default branch
-    if: github.base_ref == github.event.repository.default_branch
-    runs-on: ubuntu-latest
-    steps:
-    - uses: julia-actions/setup-julia@latest
-      with:
-        version: '1'
-    - uses: actions/checkout@v4
-    - uses: julia-actions/julia-buildpkg@latest
-    - uses: julia-actions/julia-invalidations@v1
-      id: invs_pr
-
-    - uses: actions/checkout@v4
-      with:
-        ref: ${{ github.event.repository.default_branch }}
-    - uses: julia-actions/julia-buildpkg@latest
-    - uses: julia-actions/julia-invalidations
-      id: invs_default
-    
-    - name: Report invalidation counts
-      run: |
-        echo "Invalidations on default branch: ${{ steps.invs_default.outputs.total }} (${{ steps.invs_default.outputs.deps }} via deps)" >> $GITHUB_STEP_SUMMARY
-        echo "This branch: ${{ steps.invs_pr.outputs.total }} (${{ steps.invs_pr.outputs.deps }} via deps)" >> $GITHUB_STEP_SUMMARY
-    - name: Check if the PR does increase number of invalidations
-      if: steps.invs_pr.outputs.total > steps.invs_default.outputs.total
-      run: exit 1
\ No newline at end of file
diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml
new file mode 100644
index 00000000..604011f0
--- /dev/null
+++ b/.github/workflows/SpellCheck.yml
@@ -0,0 +1,13 @@
+name: Spell Check
+
+on: [pull_request]
+
+jobs:
+  typos-check:
+    name: Spell Check with Typos
+    runs-on: ubuntu-latest
+    steps:
+      - name: Checkout Actions Repository
+        uses: actions/checkout@v4
+      - name: Check spelling
+        uses: crate-ci/typos@v1.27.0
\ No newline at end of file
diff --git a/.typos.toml b/.typos.toml
new file mode 100644
index 00000000..6fb5b8e5
--- /dev/null
+++ b/.typos.toml
@@ -0,0 +1,3 @@
+[default.extend-words]
+# Don't correct the surname "Teh"
+multline = "multline"
\ No newline at end of file
diff --git a/docs/make.jl b/docs/make.jl
index 5ea4d394..08a90f6a 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -5,7 +5,7 @@ using Documenter
 using DocumenterVitepress
 using DocumenterCitations
 
-# extentions
+# extensions
 using ModelingToolkit
 using OrdinaryDiffEqTsit5
 using SteadyStateDiffEq
diff --git a/docs/src/background/limit_cycles.md b/docs/src/background/limit_cycles.md
index fbc63f1e..bba743fa 100644
--- a/docs/src/background/limit_cycles.md
+++ b/docs/src/background/limit_cycles.md
@@ -1,6 +1,6 @@
 # [Limit cycles](@id limit_cycles_bg)
 
-We explain how HarmonicBalance.jl uses a new technique to find limit cycles in systems of nonlinear ODEs. For a more in depth overwiew see Chapter 6 in [Jan Košata's PhD theses](https://www.doi.org/10.3929/ethz-b-000589190) or [del_Pino_2024](https://www.doi.org/10.1103/PhysRevResearch.6.03318). 
+We explain how HarmonicBalance.jl uses a new technique to find limit cycles in systems of nonlinear ODEs. For a more in depth overview see Chapter 6 in [Jan Košata's PhD theses](https://www.doi.org/10.3929/ethz-b-000589190) or [del_Pino_2024](https://www.doi.org/10.1103/PhysRevResearch.6.03318). 
 
 ## Limit cycles from a Hopf bifurcation
 The end product of the [harmonic balance technique](@ref intro_hb) are what we call the harmonic equations, i.e., first-order ODEs for the harmonic variables $\mathbf{U}(T)$:
diff --git a/docs/src/examples/parametron.md b/docs/src/examples/parametron.md
index 4db5b685..aef1ff02 100644
--- a/docs/src/examples/parametron.md
+++ b/docs/src/examples/parametron.md
@@ -52,7 +52,7 @@ and replacing this by the time independent (averaged) equations of motion. This
 harmonic_eq = get_harmonic_equations(diff_eq)
 ````
 
-The output of these equations are consistent with the result found in the literature. Now we are interested in the linear response spectrum, which we can obtain from the solutions to the averaged equations (rotating frame) as a function of the external drive, after fixing all other parameters in the system. A call to `get_steady_states` then retrieves all steadystates found allong the sweep employing the homotopy continuation method, which occurs in a complex space (see the nice [HomotopyContinuation.jl docs](https://www.juliahomotopycontinuation.org))
+The output of these equations are consistent with the result found in the literature. Now we are interested in the linear response spectrum, which we can obtain from the solutions to the averaged equations (rotating frame) as a function of the external drive, after fixing all other parameters in the system. A call to `get_steady_states` then retrieves all steadystates found along the sweep employing the homotopy continuation method, which occurs in a complex space (see the nice [HomotopyContinuation.jl docs](https://www.juliahomotopycontinuation.org))
 
 ## 1D parameters
 We start with a `varied` set containing one parameter, $\omega$,
@@ -64,7 +64,7 @@ varied = ω => range(0.9, 1.1, 100)
 result = get_steady_states(harmonic_eq, varied, fixed)
 ````
 
-In `get_steady_states`, the default method `WarmUp()` initiates the homotopy in a generalised version of the harmonic equations, where parameters become random complex numbers. A parameter homotopy then follows to each of the frequency values $\omega$ in sweep. This offers speed-up, but requires to be tested in each scenario againts the method `TotalDegree`, which initializes the homotopy in a total degree system (maximum number of roots), but needs to track significantly more homotopy paths and there is slower.
+In `get_steady_states`, the default method `WarmUp()` initiates the homotopy in a generalised version of the harmonic equations, where parameters become random complex numbers. A parameter homotopy then follows to each of the frequency values $\omega$ in sweep. This offers speed-up, but requires to be tested in each scenario against the method `TotalDegree`, which initializes the homotopy in a total degree system (maximum number of roots), but needs to track significantly more homotopy paths and there is slower.
 
 After solving the system, we can save the full output of the simulation and the model (e.g. symbolic expressions for the harmonic equations) into a file
 
@@ -73,13 +73,13 @@ HarmonicBalance.save("parametron_result.jld2", result);
 nothing #hide
 ````
 
-During the execution of `get_steady_states`, different solution branches are classified by their proximity in complex space, with subsequent filtering of real (physically accceptable solutions). In addition, the stability properties of each steady state is assesed from the eigenvalues of the Jacobian matrix. All this information can be succintly represented in a 1D plot via
+During the execution of `get_steady_states`, different solution branches are classified by their proximity in complex space, with subsequent filtering of real (physically acceptable solutions). In addition, the stability properties of each steady state is assessed from the eigenvalues of the Jacobian matrix. All this information can be succinctly represented in a 1D plot via
 
 ````@example parametron
 plot(result; x="ω", y="sqrt(u1^2 + v1^2)")
 ````
 
-The user can also introduce custom clases based on parameter conditions via `classify_solutions!`. Plots can be overlaid and use keywords from `Plots`,\
+The user can also introduce custom classes based on parameter conditions via `classify_solutions!`. Plots can be overlaid and use keywords from `Plots`,\
 
 ````@example parametron
 classify_solutions!(result, "sqrt(u1^2 + v1^2) > 0.1", "large")
diff --git a/docs/src/examples/wave_mixing.md b/docs/src/examples/wave_mixing.md
index 2be86998..0c61b42e 100644
--- a/docs/src/examples/wave_mixing.md
+++ b/docs/src/examples/wave_mixing.md
@@ -52,7 +52,7 @@ plot(p1, p2, p3; layout=(1, 3), size=(900, 300), margin=5mm)
 If we only have a cubic nonlineariy $\alpha$, we observe the normal duffing oscillator
 response with no response at $2\omega$.
 
-We would like to investigate the three-wave mixing of the driven Duffing oscillator. This means we can excite the system resonantly if the oscillation frequencies $\omega_1$ and $\omega_2$ fullfil the conditions $\omega_1\pm\omega_2=\pm\omega_0$. Here, we will especially focus on the degenerate three wave mixing, where $\omega_2=\omega_0$ such that $2\omega_0=\omega_1$. This is a very important process in quantum optics, since it allows us to generate photons with a frequency in the visible range from photons with a frequency in the infrared range. This is called frequency doubling and is used in many applications, e.g. in laser pointers.
+We would like to investigate the three-wave mixing of the driven Duffing oscillator. This means we can excite the system resonantly if the oscillation frequencies $\omega_1$ and $\omega_2$ fulfil the conditions $\omega_1\pm\omega_2=\pm\omega_0$. Here, we will especially focus on the degenerate three wave mixing, where $\omega_2=\omega_0$ such that $2\omega_0=\omega_1$. This is a very important process in quantum optics, since it allows us to generate photons with a frequency in the visible range from photons with a frequency in the infrared range. This is called frequency doubling and is used in many applications, e.g. in laser pointers.
 
 ````@example wave_mixing
 varied = (ω => range(0.9, 1.2, 200))
@@ -70,7 +70,7 @@ plot(p1, p2, p3; layout=(1, 3), size=(900, 300), margin=5mm)
 If we only have a cubic nonlineariy $\alpha$, we observe the normal duffing oscillator
 response with no response at $2\omega$.
 
-We would like to investigate the three-wave mixing of the driven Duffing oscillator. This means we can excite the system resonantly if the oscillation frequencies $\omega_1$ and $\omega_2$ fullfil the conditions $\omega_1\pm\omega_2=\pm\omega_0$. Here, we will especially focus on the degenerate three wave mixing, where $\omega_2=\omega_0$ such that $2\omega_0=\omega_1$. This is a very important process in quantum optics, since it allows us to generate photons with a frequency in the visible range from photons with a frequency in the infrared range. This is called frequency doubling and is used in many applications, e.g. in laser pointers.
+We would like to investigate the three-wave mixing of the driven Duffing oscillator. This means we can excite the system resonantly if the oscillation frequencies $\omega_1$ and $\omega_2$ fulfil the conditions $\omega_1\pm\omega_2=\pm\omega_0$. Here, we will especially focus on the degenerate three wave mixing, where $\omega_2=\omega_0$ such that $2\omega_0=\omega_1$. This is a very important process in quantum optics, since it allows us to generate photons with a frequency in the visible range from photons with a frequency in the infrared range. This is called frequency doubling and is used in many applications, e.g. in laser pointers.
 
 ````@example wave_mixing
 varied = (ω => range(0.9, 1.2, 200))
diff --git a/docs/src/manual/extracting_harmonics.md b/docs/src/manual/extracting_harmonics.md
index 8e502b36..653ff84d 100644
--- a/docs/src/manual/extracting_harmonics.md
+++ b/docs/src/manual/extracting_harmonics.md
@@ -4,7 +4,7 @@
 Once a `DifferentialEquation` is defined and its harmonics specified, one can extract the harmonic equations using `get_harmonic_equations`, which itself is composed of the subroutines `harmonic_ansatz`, `slow_flow`, `fourier_transform!` and `drop_powers`. 
 
 The harmonic equations use an additional time variable specified as `slow_time` in `get_harmonic_equations`. This is essentially a label distinguishing the time dependence of the harmonic variables (expected to be slow)
-from that of the oscillating terms (expeted to be fast). When the equations are Fourier-transformed to remove oscillating terms, `slow_time` is treated as a constant. Such an approach is exact when looking for steady states. 
+from that of the oscillating terms (expected to be fast). When the equations are Fourier-transformed to remove oscillating terms, `slow_time` is treated as a constant. Such an approach is exact when looking for steady states. 
 
 ```@docs
 get_harmonic_equations
diff --git a/docs/src/manual/time_dependent.md b/docs/src/manual/time_dependent.md
index e86e8632..01942912 100644
--- a/docs/src/manual/time_dependent.md
+++ b/docs/src/manual/time_dependent.md
@@ -2,7 +2,7 @@
 
 Generally, solving the ODE of oscillatory systems in time requires numerically tracking the oscillations. This is a computationally expensive process; however, using the harmonic ansatz removes the oscillatory time-dependence. Simulating instead the harmonic variables of a `HarmonicEquation` is vastly more efficient - a steady state of the system appears as a fixed point in multidimensional space rather than an oscillatory function.
 
-The Extention `TimeEvolution` is used to interface `HarmonicEquation` with the solvers contained in `OrdinaryDiffEq.jl`. Time-dependent parameter sweeps are defined using the object `AdiabaticSweep`. To use the `TimeEvolution` extension, one must first load the `OrdinaryDiffEq.jl` package.
+The extension `TimeEvolution` is used to interface `HarmonicEquation` with the solvers contained in `OrdinaryDiffEq.jl`. Time-dependent parameter sweeps are defined using the object `AdiabaticSweep`. To use the `TimeEvolution` extension, one must first load the `OrdinaryDiffEq.jl` package.
 ```@docs
 ODEProblem(::HarmonicEquation, ::Any; timespan::Tuple)
 AdiabaticSweep
diff --git a/docs/src/tutorials/classification.md b/docs/src/tutorials/classification.md
index b79178ec..c384060a 100644
--- a/docs/src/tutorials/classification.md
+++ b/docs/src/tutorials/classification.md
@@ -15,7 +15,7 @@ add_harmonic!(diff_eq, x, ω);
 harmonic_eq = get_harmonic_equations(diff_eq)
 ```
 
-We performe a 2d sweep in the driving frequency $\omega$ and driving strength $\lambda$:
+We perform a 2d sweep in the driving frequency $\omega$ and driving strength $\lambda$:
 ```@example classification
 fixed = (ω₀ => 1.0, γ => 0.002, α => 1.0)
 varied = (ω => range(0.99, 1.01, 100), λ => range(1e-6, 0.03, 100))
@@ -39,7 +39,7 @@ plot(result_2D, y="√(u1^2+v1^2)", cut=λ => 0.01, class="stable") |> display
 ```@example classification
 get_single_solution(result_2D; branch=1, index=(1, 1))
 ```
-This solution becomes stable again outside the green lobe. Also called Mathieu lobe. Indeed, we can classify the zero amplitude solution by adding an extra catagory as a class:
+This solution becomes stable again outside the green lobe. Also called Mathieu lobe. Indeed, we can classify the zero amplitude solution by adding an extra category as a class:
 ```@example classification
 classify_solutions!(result_2D, "sqrt(u1^2 + v1^2) < 0.001", "zero")
 result_2D
diff --git a/docs/src/tutorials/limit_cycles.md b/docs/src/tutorials/limit_cycles.md
index acd790a9..aa3face3 100644
--- a/docs/src/tutorials/limit_cycles.md
+++ b/docs/src/tutorials/limit_cycles.md
@@ -47,7 +47,7 @@ So far, we have largely focused on finding and analysing steady states, i.e., fi
 \end{equation}
 ```
 
-Fixed points are however merely a subset of possible solutions of Eq. \eqref{eq:harmeqfull} -- strictly speaking, solutions where $\mathbf{u}(T)$ remains time-dependent are allowed. These are quite unusual, since $\bar{\mathbf{F}} (\mathbf{u})$ [is by construction time-independent](@ref intro_hb) and Eq. \eqref{eq:harmeqfull} thus possesses _continuous time-translation symmetry_. The appearance of explicitly time-dependent solutions then consitutes spontaneous time-translation symmetry breaking.
+Fixed points are however merely a subset of possible solutions of Eq. \eqref{eq:harmeqfull} -- strictly speaking, solutions where $\mathbf{u}(T)$ remains time-dependent are allowed. These are quite unusual, since $\bar{\mathbf{F}} (\mathbf{u})$ [is by construction time-independent](@ref intro_hb) and Eq. \eqref{eq:harmeqfull} thus possesses _continuous time-translation symmetry_. The appearance of explicitly time-dependent solutions then constitutes spontaneous time-translation symmetry breaking.
 
 Such solutions, known as _limit cycles_, typically appear as closed periodic trajectories of the harmonic variables $\mathbf{u}(T)$. The simplest way to numerically characterise them is a time-dependent simulation, using a steady-state diagram as a guide.
 
diff --git a/docs/src/tutorials/linear_response.md b/docs/src/tutorials/linear_response.md
index 1c3ab2de..79ad52a5 100644
--- a/docs/src/tutorials/linear_response.md
+++ b/docs/src/tutorials/linear_response.md
@@ -78,7 +78,7 @@ plot(
     size=(600, 250), margin=3mm
 )
 ```
-In branch 1 the linear response to white noise shows _more than one peak_. This is a distinctly nonlinear phenomenon, indicitive if the squeezing of the steady state. Branch 2 is again quasi-linear, which stems from its low amplitude.
+In branch 1 the linear response to white noise shows _more than one peak_. This is a distinctly nonlinear phenomenon, indicative if the squeezing of the steady state. Branch 2 is again quasi-linear, which stems from its low amplitude.
 
 Following [Huber et al.](https://doi.org/10.1103/PhysRevX.10.021066), we may also fix $\omega = \omega_0$ and plot the linear response as a function of $F$. The response turns out to be single-valued over a large range of driving strengths. Using a log scale for the x-axis:
 
diff --git a/docs/src/tutorials/steady_states.md b/docs/src/tutorials/steady_states.md
index 0eea3b76..797ce8e7 100644
--- a/docs/src/tutorials/steady_states.md
+++ b/docs/src/tutorials/steady_states.md
@@ -45,7 +45,7 @@ Now everything is ready to crank the handle. `get_steady_states` solves our `har
 ```@example steady_state
 result = get_steady_states(harmonic_eq, varied, fixed)
 ```
-The algorithm has found 3 solution branches in total (out of the [hypothetically admissible](https://en.wikipedia.org/wiki/B%C3%A9zout%27s_theorem) ``3^{2} = 9``). All of these are real -- and thefore physically observable -- for at least some values of ``\omega``. Only 2 branches are stable under infinitesimal perturbations. The "Classes" are boolean labels classifying each solution point, which may be used to select results for plotting.
+The algorithm has found 3 solution branches in total (out of the [hypothetically admissible](https://en.wikipedia.org/wiki/B%C3%A9zout%27s_theorem) ``3^{2} = 9``). All of these are real -- and therefore physically observable -- for at least some values of ``\omega``. Only 2 branches are stable under infinitesimal perturbations. The "Classes" are boolean labels classifying each solution point, which may be used to select results for plotting.
 
 We now want to visualize the results. Here we plot the solution amplitude, ``\sqrt{U^2 + V^2}`` against the drive frequency ``\omega``: 
 ```@example steady_state
diff --git a/examples/parametron.jl b/examples/parametron.jl
index 0db5a953..6636e3f7 100644
--- a/examples/parametron.jl
+++ b/examples/parametron.jl
@@ -40,7 +40,7 @@ nothing # hide
 
 harmonic_eq = get_harmonic_equations(diff_eq)
 
-# The output of these equations are consistent with the result found in the literature. Now we are interested in the linear response spectrum, which we can obtain from the solutions to the averaged equations (rotating frame) as a function of the external drive, after fixing all other parameters in the system. A call to `get_steady_states` then retrieves all steadystates found allong the sweep employing the homotopy continuation method, which occurs in a complex space (see the nice [HomotopyContinuation.jl docs](https://www.juliahomotopycontinuation.org))
+# The output of these equations are consistent with the result found in the literature. Now we are interested in the linear response spectrum, which we can obtain from the solutions to the averaged equations (rotating frame) as a function of the external drive, after fixing all other parameters in the system. A call to `get_steady_states` then retrieves all steadystates found along the sweep employing the homotopy continuation method, which occurs in a complex space (see the nice [HomotopyContinuation.jl docs](https://www.juliahomotopycontinuation.org))
 
 # ## 1D parameters
 # We start with a `varied` set containing one parameter, $\omega$,
@@ -50,17 +50,17 @@ varied = ω => range(0.9, 1.1, 100)
 
 result = get_steady_states(harmonic_eq, varied, fixed)
 
-# In `get_steady_states`, the default method `WarmUp()` initiates the homotopy in a generalised version of the harmonic equations, where parameters become random complex numbers. A parameter homotopy then follows to each of the frequency values $\omega$ in sweep. This offers speed-up, but requires to be tested in each scenario againts the method `TotalDegree`, which initializes the homotopy in a total degree system (maximum number of roots), but needs to track significantly more homotopy paths and there is slower.
+# In `get_steady_states`, the default method `WarmUp()` initiates the homotopy in a generalised version of the harmonic equations, where parameters become random complex numbers. A parameter homotopy then follows to each of the frequency values $\omega$ in sweep. This offers speed-up, but requires to be tested in each scenario against the method `TotalDegree`, which initializes the homotopy in a total degree system (maximum number of roots), but needs to track significantly more homotopy paths and there is slower.
 
 # After solving the system, we can save the full output of the simulation and the model (e.g. symbolic expressions for the harmonic equations) into a file
 
 HarmonicBalance.save("parametron_result.jld2", result);
 
-# During the execution of `get_steady_states`, different solution branches are classified by their proximity in complex space, with subsequent filtering of real (physically accceptable solutions). In addition, the stability properties of each steady state is assesed from the eigenvalues of the Jacobian matrix. All this information can be succintly represented in a 1D plot via
+# During the execution of `get_steady_states`, different solution branches are classified by their proximity in complex space, with subsequent filtering of real (physically acceptable solutions). In addition, the stability properties of each steady state is assessed from the eigenvalues of the Jacobian matrix. All this information can be succinctly represented in a 1D plot via
 
 plot(result; x="ω", y="sqrt(u1^2 + v1^2)")
 
-# The user can also introduce custom clases based on parameter conditions via `classify_solutions!`. Plots can be overlaid and use keywords from `Plots`,\
+# The user can also introduce custom classes based on parameter conditions via `classify_solutions!`. Plots can be overlaid and use keywords from `Plots`,\
 
 classify_solutions!(result, "sqrt(u1^2 + v1^2) > 0.1", "large")
 plot(result, "sqrt(u1^2 + v1^2)"; class=["physical", "large"], style=:dash)
diff --git a/examples/wave_mixing.jl b/examples/wave_mixing.jl
index db548ac0..a8a50505 100644
--- a/examples/wave_mixing.jl
+++ b/examples/wave_mixing.jl
@@ -42,7 +42,7 @@ plot(p1, p2, p3; layout=(1, 3), size=(900, 300), margin=5mm)
 # If we only have a cubic nonlineariy $\alpha$, we observe the normal duffing oscillator
 # response with no response at $2\omega$.
 
-# We would like to investigate the three-wave mixing of the driven Duffing oscillator. This means we can excite the system resonantly if the oscillation frequencies $\omega_1$ and $\omega_2$ fullfil the conditions $\omega_1\pm\omega_2=\pm\omega_0$. Here, we will especially focus on the degenerate three wave mixing, where $\omega_2=\omega_0$ such that $2\omega_0=\omega_1$. This is a very important process in quantum optics, since it allows us to generate photons with a frequency in the visible range from photons with a frequency in the infrared range. This is called frequency doubling and is used in many applications, e.g. in laser pointers.
+# We would like to investigate the three-wave mixing of the driven Duffing oscillator. This means we can excite the system resonantly if the oscillation frequencies $\omega_1$ and $\omega_2$ fulfil the conditions $\omega_1\pm\omega_2=\pm\omega_0$. Here, we will especially focus on the degenerate three wave mixing, where $\omega_2=\omega_0$ such that $2\omega_0=\omega_1$. This is a very important process in quantum optics, since it allows us to generate photons with a frequency in the visible range from photons with a frequency in the infrared range. This is called frequency doubling and is used in many applications, e.g. in laser pointers.
 
 varied = (ω => range(0.9, 1.2, 200))
 fixed = (α => 0.0, β => 1.0, ω0 => 1.0, γ => 0.005, F => 0.0025)
@@ -58,7 +58,7 @@ plot(p1, p2, p3; layout=(1, 3), size=(900, 300), margin=5mm)
 # If we only have a cubic nonlineariy $\alpha$, we observe the normal duffing oscillator
 # response with no response at $2\omega$.
 
-# We would like to investigate the three-wave mixing of the driven Duffing oscillator. This means we can excite the system resonantly if the oscillation frequencies $\omega_1$ and $\omega_2$ fullfil the conditions $\omega_1\pm\omega_2=\pm\omega_0$. Here, we will especially focus on the degenerate three wave mixing, where $\omega_2=\omega_0$ such that $2\omega_0=\omega_1$. This is a very important process in quantum optics, since it allows us to generate photons with a frequency in the visible range from photons with a frequency in the infrared range. This is called frequency doubling and is used in many applications, e.g. in laser pointers.
+# We would like to investigate the three-wave mixing of the driven Duffing oscillator. This means we can excite the system resonantly if the oscillation frequencies $\omega_1$ and $\omega_2$ fulfil the conditions $\omega_1\pm\omega_2=\pm\omega_0$. Here, we will especially focus on the degenerate three wave mixing, where $\omega_2=\omega_0$ such that $2\omega_0=\omega_1$. This is a very important process in quantum optics, since it allows us to generate photons with a frequency in the visible range from photons with a frequency in the infrared range. This is called frequency doubling and is used in many applications, e.g. in laser pointers.
 
 varied = (ω => range(0.9, 1.2, 200))
 fixed = (α => 1.0, β => 1.0, ω0 => 1.0, γ => 0.005, F => 0.0025)
diff --git a/ext/ModelingToolkitExt.jl b/ext/ModelingToolkitExt.jl
index ccfbabaa..63afe8f9 100644
--- a/ext/ModelingToolkitExt.jl
+++ b/ext/ModelingToolkitExt.jl
@@ -108,7 +108,7 @@ function ModelingToolkit.NonlinearProblem(
     eom::HarmonicEquation, u0, p::ParameterList; in_place=true, kwargs...
 )
     ss_prob = SteadyStateProblem(eom, u0, p::ParameterList; in_place=in_place, kwargs...)
-    return NonlinearProblem(ss_prob) # gives warning of some internal deprication
+    return NonlinearProblem(ss_prob)
 end
 
 function ModelingToolkit.SteadyStateProblem(
diff --git a/src/ExprUtils/drop_powers.jl b/src/ExprUtils/drop_powers.jl
index 7c28b101..1619c794 100644
--- a/src/ExprUtils/drop_powers.jl
+++ b/src/ExprUtils/drop_powers.jl
@@ -40,7 +40,7 @@ end
 drop_powers(expr, var::Num, deg::Int) = drop_powers(expr, [var], deg)
 drop_powers(x, vars, deg::Int) = drop_powers(Num(x), vars, deg)
 
-"Return the highest power of `y` occuring in the term `x`."
+"Return the highest power of `y` occurring in the term `x`."
 function max_power(x::Num, y::Num)
     terms = get_all_terms(x)
     powers = power_of.(terms, y)
diff --git a/src/FFTWExt.jl b/src/FFTWExt.jl
index 370cc4c2..2ca57dda 100644
--- a/src/FFTWExt.jl
+++ b/src/FFTWExt.jl
@@ -5,7 +5,7 @@ using FFTW: fft, fftfreq, fftshift
 using Peaks: Peaks
 
 """
-Fourier transform the timeseries of a simulation in the rotating frame and calculate the quadratures and freqeuncies in the non-rotating frame.
+Fourier transform the timeseries of a simulation in the rotating frame and calculate the quadratures and frequencies in the non-rotating frame.
 """
 function FFT(soln_u, soln_t; window=DSP.Windows.hanning)
     "Input: solution object of DifferentialEquation (positions array and corresponding time)
diff --git a/src/HC_wrapper/homotopy_interface.jl b/src/HC_wrapper/homotopy_interface.jl
index cde9bc61..83b37c8b 100644
--- a/src/HC_wrapper/homotopy_interface.jl
+++ b/src/HC_wrapper/homotopy_interface.jl
@@ -61,7 +61,7 @@ function HarmonicBalance.Problem(
     system = HomotopyContinuation.System(eqs_HC; variables=conv_vars, parameters=conv_para)
     J = HarmonicBalance.get_Jacobian(equations, variables) #all derivatives are assumed to be in the left hand side;
     return Problem(variables, parameters, system, J)
-end # TODO is this funciton still needed/used?
+end # TODO is this function still needed/used?
 
 function System(eom::HarmonicEquation)
     eqs = expand_derivatives.(_remove_brackets(eom))
diff --git a/src/HarmonicBalance.jl b/src/HarmonicBalance.jl
index f8d22f61..4ffedeb2 100644
--- a/src/HarmonicBalance.jl
+++ b/src/HarmonicBalance.jl
@@ -38,7 +38,7 @@ using SymbolicUtils: SymbolicUtils
 include("ExprUtils/ExprUtils.jl")
 using .ExprUtils: is_harmonic, substitute_all, drop_powers, count_derivatives
 
-include("extention_functions.jl")
+include("extension_functions.jl")
 include("utils.jl")
 include("types.jl")
 include("DifferentialEquation.jl")
diff --git a/src/HarmonicEquation.jl b/src/HarmonicEquation.jl
index 1a126851..79020cb9 100644
--- a/src/HarmonicEquation.jl
+++ b/src/HarmonicEquation.jl
@@ -72,7 +72,7 @@ function harmonic_ansatz(diff_eom::DifferentialEquation, time::Num)
     a_idx = 1
 
     for nvar in get_variables(diff_eom) # sum over natural variables
-        to_substitute = Num(0) # combine all the subtitution rules for var
+        to_substitute = Num(0) # combine all the substitution rules for var
         for ω in diff_eom.harmonics[nvar]
             if !isequal(ω, 0) # nonzero harmonic - create u,v
                 rule_u, hvar_u = _create_harmonic_variable(
diff --git a/src/KrylovBogoliubov/KrylovEquation.jl b/src/KrylovBogoliubov/KrylovEquation.jl
index eb308be0..8f8eb35a 100644
--- a/src/KrylovBogoliubov/KrylovEquation.jl
+++ b/src/KrylovBogoliubov/KrylovEquation.jl
@@ -45,7 +45,7 @@ function get_krylov_equations(
     diff_eom::DifferentialEquation; order, fast_time=nothing, slow_time=nothing
 )
     order < 1 && error("The order of the Krylov-Bogoliubov method must be at least 1!")
-    order > 2 && error("Krylov-Bogoliubov implemetation only supports up to second order!")
+    order > 2 && error("Krylov-Bogoliubov implementation only supports up to second order!")
 
     slow_time = isnothing(slow_time) ? (@variables T; T) : slow_time
     fast_time = isnothing(fast_time) ? get_independent_variables(diff_eom)[1] : fast_time
@@ -53,14 +53,14 @@ function get_krylov_equations(
     harmonics = values(diff_eom.harmonics)
     all(isempty.(harmonics)) && error("No harmonics specified!")
     any(isempty.(harmonics)) &&
-        error("Krylov-Bogoliubov method needs all vairables to have a single harmonic!")
+        error("Krylov-Bogoliubov method needs all variables to have a single harmonic!")
     any(length.(harmonics) .> 1) &&
         error("Krylov-Bogoliubov method only supports a single harmonic!")
 
-    deom = deepcopy(diff_eom)
-    !is_rearranged_standard(deom) ? rearrange_standard!(deom) : nothing
-    first_order_transform!(deom, fast_time)
-    eom = van_der_Pol(deom, fast_time)
+    dEOM = deepcopy(diff_eom)
+    !is_rearranged_standard(dEOM) ? rearrange_standard!(dEOM) : nothing
+    first_order_transform!(dEOM, fast_time)
+    eom = van_der_Pol(dEOM, fast_time)
 
     eom = slow_flow(eom; fast_time=fast_time, slow_time=slow_time, degree=2)
 
diff --git a/src/LinearResponse/jacobians.jl b/src/LinearResponse/jacobians.jl
index 706d8158..ff753c91 100644
--- a/src/LinearResponse/jacobians.jl
+++ b/src/LinearResponse/jacobians.jl
@@ -45,7 +45,7 @@ function get_Jacobian(eqs::Vector{Equation}, vars::Vector{Num})
 end
 
 """
-Code folllows for an implicit treatment of the Jacobian. Usually we rearrange the linear response equations to have time-derivatives on one side.
+Code follows for an implicit treatment of the Jacobian. Usually we rearrange the linear response equations to have time-derivatives on one side.
     This may be extremely costly. Implicit evaluation means only solving the equations AFTER numerical values have been plugged in, giving
     a constant time cost per run.
 """
diff --git a/src/LinearResponse/plotting.jl b/src/LinearResponse/plotting.jl
index a947aa90..d8d31ca9 100644
--- a/src/LinearResponse/plotting.jl
+++ b/src/LinearResponse/plotting.jl
@@ -214,7 +214,7 @@ end
 """
     plot_rotframe_jacobian_response(res::Result; Ω_range, branch::Int, logscale=true, damping_mod::Float64 = 1.0, show_progress=true, kwargs...)
 
-Plot the linear response to white noise in the rotating frame for Result `res` on `branch` for input frequencies `Ω_range`. 'damping_mod' gets multiplied by the real part of the eigenvalues of the Jacobian in order to be able to make peaks with similar frequency seperately identifiable.
+Plot the linear response to white noise in the rotating frame for Result `res` on `branch` for input frequencies `Ω_range`. 'damping_mod' gets multiplied by the real part of the eigenvalues of the Jacobian in order to be able to make peaks with similar frequency separately identifiable.
 
 Any kwargs are fed to Plots' gr().
 
diff --git a/src/extention_functions.jl b/src/extension_functions.jl
similarity index 100%
rename from src/extention_functions.jl
rename to src/extension_functions.jl
diff --git a/src/plotting_Plots.jl b/src/plotting_Plots.jl
index 3f5bc161..71b04c2e 100644
--- a/src/plotting_Plots.jl
+++ b/src/plotting_Plots.jl
@@ -218,7 +218,7 @@ function plot2D_cut(
     end
 
     cut_par, cut_value = cut
-    # compare strings beacuse type Num cannot be compared
+    # compare strings because type Num cannot be compared
     swept_pars = res.swept_parameters.keys
     x_index = findfirst(sym -> string(sym) != string(cut_par), swept_pars)
     isnothing(x_index) && error("The variable $cut_par was not swept over.")
diff --git a/src/solve_homotopy.jl b/src/solve_homotopy.jl
index f5e198ff..2815d40a 100644
--- a/src/solve_homotopy.jl
+++ b/src/solve_homotopy.jl
@@ -36,7 +36,7 @@ A steady state result for 100 parameter points
 
 ```
 
-It is also possible to perfrom 2-dimensional sweeps.
+It is also possible to perform 2-dimensional sweeps.
 ```julia-repl
 # The swept parameters take precedence over fixed -> use the same fixed
 julia> range = ParameterRange(ω => range(0.8,1.2,100), F => range(0.1,1.0,10) )
diff --git a/src/sorting.jl b/src/sorting.jl
index f3ddeb30..a104c90a 100644
--- a/src/sorting.jl
+++ b/src/sorting.jl
@@ -33,7 +33,7 @@ end
 ###
 
 """
-Removes rows and columns with given idices from the Matrix M.
+Removes rows and columns with given indices from the Matrix M.
 The row/column indices are defined with respect to the original M!
 """
 function remove_rows_columns(M::Matrix, rows::Vector{Int64}, cols::Vector{Int64})
diff --git a/test/runtests.jl b/test/runtests.jl
index 2a835ac9..302cfd29 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -80,15 +80,15 @@ end
     include("limit_cycle.jl")
 end
 
-@testset "Extentions" begin
-    @testset "Time evolution extention" begin
+@testset "extensions" begin
+    @testset "Time evolution extension" begin
         include("time_evolution.jl")
         include("hysteresis_sweep.jl")
     end
-    @testset "ModelingToolkit extention" begin
+    @testset "ModelingToolkit extension" begin
         include("ModelingToolkitExt.jl")
     end
-    @testset "SteadyState Extention" begin
+    @testset "SteadyState extension" begin
         include("SteadyStateDiffEqExt.jl")
     end
 end