build: add spell checker

.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/[email protected]

.typos.toml

@@ -0,0 +1,3 @@
+[default.extend-words]
+# Don't correct the surname "Teh"
+multline = "multline"

docs/make.jl

@@ -5,7 +5,7 @@ using Documenter
 using DocumenterVitepress
 using DocumenterCitations
 
-# extentions
+# extensions
 using ModelingToolkit
 using OrdinaryDiffEqTsit5
 using SteadyStateDiffEq

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)$:

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")

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))

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

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

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

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.
 

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:
 

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