diff --git a/README.md b/README.md
index 3fa6b76e..6a777451 100644
--- a/README.md
+++ b/README.md
@@ -10,7 +10,7 @@
 [![JET](https://img.shields.io/badge/%E2%9C%88%EF%B8%8F%20tested%20with%20-%20JET.jl%20-%20red)](https://github.com/aviatesk/JET.jl)
 
 
-**HarmonicBalance.jl** is a Julia package for solving nonlinear differential equations using the method of harmonic balance.
+**HarmonicBalance.jl** is a Julia package for solving periodic differential equations using the method of harmonic balance.
 
 ## Installation
 
diff --git a/docs/src/.vitepress/config.mts b/docs/src/.vitepress/config.mts
index b2473222..268c9eaf 100644
--- a/docs/src/.vitepress/config.mts
+++ b/docs/src/.vitepress/config.mts
@@ -76,8 +76,8 @@ export default defineConfig({
     sidebar: {
       "/introduction/": [
         { text: 'Introduction', link: '/index' },
-        { text: 'Overview', link: '/introduction/overview' },
-        { text: 'Tutorials', link: '/tutorials/index' },
+        // { text: 'Overview', link: '/introduction/overview' },
+        // { text: 'Tutorials', link: '/tutorials/index' },
         { text: 'Citation', link: '/introduction/citation' }
       ],
       "/background/": [
diff --git a/docs/src/background/stability_response.md b/docs/src/background/stability_response.md
index d3eefd43..37278539 100644
--- a/docs/src/background/stability_response.md
+++ b/docs/src/background/stability_response.md
@@ -2,7 +2,7 @@
 
 The core of the harmonic balance method is expressing the system's behaviour in terms of Fourier components or _harmonics_. For an $N$-coordinate system, we _choose_ a set of $M_i$ harmonics to describe each coordinate $x_i$ : 
 ```math
-\begin{equation} \label{eq:ansatz}
+\begin{equation}
 x_i(t) = \sum_{j=1}^{M_i} u_{i,j}  (T)  \cos(\omega_{i,j} t)+ v_{i,j} (T) \sin(\omega_{i,j} t) \:,
 \end{equation}
 ```
@@ -24,13 +24,13 @@ where $\bar{\mathbf{F}}(\mathbf{u})$ is a nonlinear function. A steady state $\m
 
 Let us assume that we found a steady state $\mathbf{u}_0$. When the system is in this state, it responds to small perturbations either by returning to $\mathbf{u}_0$ over some characteristic timescale (_stable state_) or by evolving away from $\mathbf{u}_0$ (_unstable state_). To analyze the stability of $\mathbf{u}_0$, we linearize the equations of motion around $\mathbf{u}_0$ for a small perturbation $\delta \mathbf{u} = \mathbf{u} - \mathbf{u}_0$ to obtain
 ```math
-\begin{equation} \label{eq:Jaceq}
+\begin{equation}
 \frac{d}{dT} \left[\delta \mathbf{u}(T)\right] =  J(\mathbf{u}_0) \delta \mathbf{u}(T) \,,
 \end{equation}
 ```
 where $J(\mathbf{u}_0)=\nabla_{\mathbf{u}}  \bar{\mathbf{F}}|_{\mathbf{u}=\mathbf{u}_0}$ is the _Jacobian matrix_ of the system evaluated at $\mathbf{u}=\mathbf{u}_0$.
 
-Eq. \eqref{eq:Jaceq} is exactly solvable for $\delta \mathbf{u}(T)$ given an initial condition $\delta \mathbf{u}(T_0)$. The solution can be expanded in terms of the complex eigenvalues $\lambda_r$ and eigenvectors $\mathbf{v}_r$ of $J(\mathbf{u}_0)$, namely
+The linearised system is exactly solvable for $\delta \mathbf{u}(T)$ given an initial condition $\delta \mathbf{u}(T_0)$. The solution can be expanded in terms of the complex eigenvalues $\lambda_r$ and eigenvectors $\mathbf{v}_r$ of $J(\mathbf{u}_0)$, namely
 
 ```math
 \begin{equation} 
@@ -51,7 +51,7 @@ The response of a stable steady state to an additional oscillatory force, caused
 \end{equation}
 ```
 
-Suppose we have found an eigenvector of $J(\mathbf{u}_0)$ such that $J(\mathbf{u}) \mathbf{v} = \lambda \mathbf{v}$. To solve Eq. \eqref{eq:Jacforced}, we insert $\delta \mathbf{u}(T) = A(\Omega)\, \mathbf{v} e^{i \Omega T}$. Projecting each side of Eq. \eqref{eq:Jacforced} onto $\mathbf{v}$ gives
+Suppose we have found an eigenvector of $J(\mathbf{u}_0)$ such that $J(\mathbf{u}) \mathbf{v} = \lambda \mathbf{v}$. To solve the linearised equations of motion, we insert $\delta \mathbf{u}(T) = A(\Omega)\, \mathbf{v} e^{i \Omega T}$. Projecting each side onto $\mathbf{v}$ gives
 
 ```math
 A(\Omega) \left( i \Omega  - \lambda \right)  = \boldsymbol{\xi} \cdot \mathbf{v} \quad \implies \quad A(\Omega) = \frac{\boldsymbol{\xi} \cdot \mathbf{v} }{-\text{Re}[\lambda] + i \left( \Omega - \text{Im}[\lambda] \right)}
@@ -59,7 +59,7 @@ A(\Omega) \left( i \Omega  - \lambda \right)  = \boldsymbol{\xi} \cdot \mathbf{v
 
 We see that each eigenvalue $\lambda$ results in a linear response that is a Lorentzian centered at $\Omega = \text{Im}[\lambda]$. Effectively, the linear response matches that of a harmonic oscillator with resonance frequency $\text{Im}[\lambda]$ and damping $\text{Re}[\lambda]$.
 
-Knowing the response of the harmonic variables $\mathbf{u}(T)$, what is the corresponding behaviour of the "natural" variables $x_i(t)$ ? To find this out, we insert the perturbation back into the harmonic ansatz \eqref{eq:ansatz}. Since we require real variables, let us use $\delta \mathbf{u}(T) = A(\Omega) \left( \mathbf{v} \, e^{i \Omega T} +   \mathbf{v}^* \, e^{-i \Omega T} \right)$. Plugging this into
+Knowing the response of the harmonic variables $\mathbf{u}(T)$, what is the corresponding behaviour of the "natural" variables $x_i(t)$? To find this out, we insert the perturbation back into the harmonic ansatz. Since we require real variables, let us use $\delta \mathbf{u}(T) = A(\Omega) \left( \mathbf{v} \, e^{i \Omega T} +   \mathbf{v}^* \, e^{-i \Omega T} \right)$. Plugging this into
 ```math
 \begin{equation}
 \delta x_{i}(t) = \sum_{j=1}^{M_i} \delta{u}_{i,j}(t) \cos(\omega_{i,j} \,t) + \delta v_{i,j} (t) \sin(\omega_{i,j} \,t) 
@@ -84,7 +84,7 @@ To make sense of this, we normalize the vector $\delta \mathbf{u}$ and use norma
 L(x)_{x_0, \gamma} = \frac{1}{\left( x - x_0 \right)^2 + \gamma^2 }
 \end{equation}
 ```
-We see that all components of $\delta x_i(t)$ [Eq. \eqref{eq:deltax}] are proportional to $L(\Omega)_{\text{Im}[\lambda], \text{Re}[\lambda]}$. The first and last two summands are Lorentzians centered at $\pm \Omega$ which oscillate at $\omega_{i,j} \pm \Omega$, respectively. From this, we can extract the linear response function in Fourier space, $\chi (\tilde{\omega})$
+We see that all components of $\delta x_i(t)$ are proportional to $L(\Omega)_{\text{Im}[\lambda], \text{Re}[\lambda]}$. The first and last two summands are Lorentzians centered at $\pm \Omega$ which oscillate at $\omega_{i,j} \pm \Omega$, respectively. From this, we can extract the linear response function in Fourier space, $\chi (\tilde{\omega})$
 
 ```math
 \begin{multline}
diff --git a/docs/src/introduction/index.md b/docs/src/introduction/index.md
index 059d4de1..2c61bdd9 100644
--- a/docs/src/introduction/index.md
+++ b/docs/src/introduction/index.md
@@ -40,9 +40,9 @@ fixed = (α => 1.0, ω0 => 1.0, F => 0.01, η => 0.1)   # fixed parameters
 varied = ω => range(0.9, 1.2, 100)           # range of parameter values
 result = get_steady_states(harmonic_eq, varied, fixed)
 ```
-The found steady states can be plotted as a function of the driving frequency:
+The obtained steady states can be plotted as a function of the driving frequency:
 ```@example getting_started
-# using Plots # hide
-# theme(:dark) # hide
 plot(result, "sqrt(u1^2 + v1^2)")
 ```
+
+If you want learn more on what you can do with HarmonicBalance.jl, check out the [tutorials](@ref tutorials). We also have collected some [examples](@ref examples) of different physical systems.
\ No newline at end of file
diff --git a/docs/src/introduction/resources.md b/docs/src/introduction/resources.md
index 11d51a92..e69de29b 100644
--- a/docs/src/introduction/resources.md
+++ b/docs/src/introduction/resources.md
@@ -1,3 +0,0 @@
-
-```@bibliography
-```
diff --git a/docs/src/tutorials/index.md b/docs/src/tutorials/index.md
index 8ce9d24c..5c5265fa 100644
--- a/docs/src/tutorials/index.md
+++ b/docs/src/tutorials/index.md
@@ -39,7 +39,7 @@ const beginner = [
 </script>
 ```
 
-# Tutorials
+# [Tutorials](@id tutorials)
 
 We show the capabilities of the package by providing a series of tutorials. Examples of other systems can be found in the [examples](@ref examples) tab.