diff --git a/docs/Project.toml b/docs/Project.toml
index 704810b6..65284347 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,7 +1,9 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
 
 [compat]
 Documenter = "0.27, 1"
+FFTW = "1.7"
 SciMLOperators = "0.2, 0.3"
diff --git a/docs/src/tutorials/fftw.md b/docs/src/tutorials/fftw.md
index eb92bc8c..69bb1499 100644
--- a/docs/src/tutorials/fftw.md
+++ b/docs/src/tutorials/fftw.md
@@ -7,7 +7,7 @@ derivative of a function.
 
 ## Copy-Paste Code
 
-```
+```@example fft
 using SciMLOperators
 using LinearAlgebra, FFTW
 
@@ -23,6 +23,12 @@ k = rfftfreq(n, 2π*n/L) |> Array
 m = length(k)
 P = plan_rfft(x)
 
+fwd(u, p, t) = P * u
+bwd(u, p, t) = P \ u
+
+fwd(du, u, p, t) = mul!(du, P, u)
+bwd(du, u, p, t) = ldiv!(du, P, u)
+
 F = FunctionOperator(fwd, x, im*k;
         T=ComplexF64,
 
@@ -50,7 +56,7 @@ in the West), a common Fourier transform library. Next, we define an equispaced
 trivial example, we already know the derivative, `du`, and write it down to later test our
 FFT wrapper.
 
-```
+```@example fft_explanation
 using SciMLOperators
 using LinearAlgebra, FFTW
 
@@ -70,10 +76,10 @@ object that can be applied to inputs that are like `x` as follows: `xhat = trans
 and `LinearAlgebra.mul!(xhat, transform, x)`.  We also get `k`, the frequency modes sampled by
 our finite grid, via the function `rfftfreq`.
 
-```
+```@example fft_explanation
 k  = rfftfreq(n, 2π*n/L) |> Array
 m  = length(k)
-tr = plan_rfft(x)
+P = plan_rfft(x)
 ```
 
 Now we are ready to define our wrapper for the FFT object. To `FunctionOperator`, we
@@ -81,7 +87,12 @@ pass the in-place forward application of the transform,
 `(du,u,p,t) -> mul!(du, transform, u)`, its inverse application,
 `(du,u,p,t) -> ldiv!(du, transform, u)`, as well as input and output prototype vectors.
 
-```
+```@example fft_explanation
+fwd(u, p, t) = P * u
+bwd(u, p, t) = P \ u
+
+fwd(du, u, p, t) = mul!(du, P, u)
+bwd(du, u, p, t) = ldiv!(du, P, u)
 F = FunctionOperator(fwd, x, im*k;
         T=ComplexF64,
 
@@ -98,15 +109,12 @@ SciMLOperators. Below, we form the derivative operator, and cache it via the fun
 `cache_operator` that requires an input prototype. We can test our derivative operator
 both in-place, and out-of-place by comparing its output to the analytical derivative.
 
-```
+```@example fft_explanation
 ik = im * DiagonalOperator(k)
 Dx = F \ ik * F
 
+Dx = cache_operator(Dx, x)
+
 @show ≈(Dx * u, du; atol=1e-8)
 @show ≈(mul!(copy(u), Dx, u), du; atol=1e-8)
 ```
-
-```
-≈(Dx * u, du; atol = 1.0e-8) = true
-≈(mul!(copy(u), Dx, u), du; atol = 1.0e-8) = true
-```