diff --git a/docs/assets/predator_prey_fig.png b/docs/assets/predator_prey_fig.png
new file mode 100644
index 0000000..808dab8
Binary files /dev/null and b/docs/assets/predator_prey_fig.png differ
diff --git a/docs/make.jl b/docs/make.jl
index 103f0b6..f753f27 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -19,7 +19,7 @@ makedocs(;
         edit_link="main",
         assets=String[],
     ),
-    pages=["Home" => "index.md", "Tutorial" =>"tutorial.md", "API" => "api.md"],
+    pages=["Home" => "index.md", "Tutorial" =>"predator_prey_tutorial.md", "API" => "api.md"],
     warnonly=true,
 )
 
diff --git a/docs/src/tutorial.md b/docs/src/predator_prey_tutorial.md
similarity index 93%
rename from docs/src/tutorial.md
rename to docs/src/predator_prey_tutorial.md
index 52a9aa9..90e1975 100644
--- a/docs/src/tutorial.md
+++ b/docs/src/predator_prey_tutorial.md
@@ -1,5 +1,5 @@
 ```@meta
-EditURL = "../../examples/predator_prey_fig.jl"
+EditURL = "../../examples/predator_prey_tutorial.jl"
 ```
 
 # Tutorial: Using LinearNoiseApproximation.jl to solve a Lotka-Volterra model
@@ -7,7 +7,7 @@ EditURL = "../../examples/predator_prey_fig.jl"
 This example demonstrates how to use `LinearNoiseApproximation.jl` package to solve a Lotka-Volterra model using the Linear Noise Approximation (LNA) and compare it to the stochastic
 trajectories obtained using the Gillespie algorithm.
 
-````@example predator_prey_fig
+````@example predator_prey_tutorial
 using LinearNoiseApproximation
 using DifferentialEquations
 using Catalyst
@@ -18,14 +18,16 @@ using Plots
 ## Define the Lotka-Volterra model
 The Lotka-Volterra model describes the dynamics of predator-prey interactions in an ecosystem. It assumes that the prey population $U$ grows at a rate represented by the parameter $\alpha$ in the absence of predators, but decreases as they are consumed by the predator population $V$ at a rate determined by the interaction strength parameter, $\beta$. The predator population decreases in size if they cannot find enough prey to consume, which is represented by the mortality rate parameter, $\delta$.
 The corresponding ODE system reads
+```math
 \begin{equation}
 	\begin{aligned}
 		\frac{\mathrm{d}U}{\mathrm{d} t} & = \alpha U - \beta U V,\\
 		\frac{\mathrm{d}V}{\mathrm{d} t} & = \beta U V  - \delta V,
 	\end{aligned}\quad t\in (0, t_{\text{end}}).
 \end{equation}
+```
 
-````@example predator_prey_fig
+````@example predator_prey_tutorial
 rn = @reaction_network begin
     @parameters α β δ
     @species U(t) V(t)
@@ -38,7 +40,7 @@ end
 ## Convert the model to a JumpSystem for the Gillespie algorithm
 Using `Catalyst.jl` we can convert `ReactionSystem` to a `JumpSystem` which can be used to simulate the stochastic trajectories using the Gillespie algorithm.
 
-````@example predator_prey_fig
+````@example predator_prey_tutorial
 jumpsys = convert(JumpSystem, rn)
 
 #define the initial conditions, parameters, and time span
@@ -101,6 +103,9 @@ end
 plt
 ````
 
+save the plot
+Plots.savefig(plt, "predator_prey_fig.png")
+
 ---
 
 *This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*
diff --git a/examples/predator_prey_fig.jl b/examples/predator_prey_tutorial.jl
similarity index 97%
rename from examples/predator_prey_fig.jl
rename to examples/predator_prey_tutorial.jl
index ab0c961..b4a70b2 100644
--- a/examples/predator_prey_fig.jl
+++ b/examples/predator_prey_tutorial.jl
@@ -13,12 +13,14 @@ using Plots
 # ## Define the Lotka-Volterra model
 # The Lotka-Volterra model describes the dynamics of predator-prey interactions in an ecosystem. It assumes that the prey population $U$ grows at a rate represented by the parameter $\alpha$ in the absence of predators, but decreases as they are consumed by the predator population $V$ at a rate determined by the interaction strength parameter, $\beta$. The predator population decreases in size if they cannot find enough prey to consume, which is represented by the mortality rate parameter, $\delta$.
 # The corresponding ODE system reads
+# ```math
 # \begin{equation}
 # 	\begin{aligned}
 # 		\frac{\mathrm{d}U}{\mathrm{d} t} & = \alpha U - \beta U V,\\
 # 		\frac{\mathrm{d}V}{\mathrm{d} t} & = \beta U V  - \delta V,
 # 	\end{aligned}\quad t\in (0, t_{\text{end}}).
 # \end{equation}
+# ```
 
 rn = @reaction_network begin
     @parameters α β δ
@@ -91,3 +93,5 @@ for (i, (mean_idx, var_idx)) in enumerate(zip(mean_idxs, var_idxs))
     )
 end
 plt
+# save the plot
+# Plots.savefig(plt, "predator_prey_fig.png")