diff --git a/docs/Project.toml b/docs/Project.toml
index 8f660d9f..b29fb873 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -4,6 +4,7 @@ CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Catlab = "134e5e36-593f-5add-ad60-77f754baafbe"
 CombinatorialSpaces = "b1c52339-7909-45ad-8b6a-6e388f7c67f2"
 ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66"
+CoordRefSystems = "b46f11dc-f210-4604-bfba-323c1ec968cb"
 Decapodes = "679ab3ea-c928-4fe6-8d59-fd451142d391"
 DiagrammaticEquations = "6f00c28b-6bed-4403-80fa-30e0dc12f317"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
diff --git a/docs/make.jl b/docs/make.jl
index 94f8e1dd..e8d7e32b 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -40,6 +40,7 @@ makedocs(
   sitename  = "Decapodes.jl",
   doctest   = false,
   checkdocs = :none,
+  draft = false;
   pagesonly = true,
   linkcheck = true,
   linkcheck_ignore = [r"agupubs\.onlinelibrary\.wiley\.com", # This gives a 403 Forbidden
@@ -49,6 +50,7 @@ makedocs(
     "Overview" => "overview/overview.md",
     "Equations" => "equations/equations.md",
     "Vortices" => "navier_stokes/ns.md",
+    "Harmonics" => "harmonics/harmonics.md",
     "Cahn-Hilliard" => "ch/cahn-hilliard.md",
     "Klausmeier" => "klausmeier/klausmeier.md",
     "CISM v2.1" => "cism/cism.md",
diff --git a/docs/src/harmonics/harmonics.md b/docs/src/harmonics/harmonics.md
new file mode 100644
index 00000000..fa124fdf
--- /dev/null
+++ b/docs/src/harmonics/harmonics.md
@@ -0,0 +1,100 @@
+# Harmonics of the Sphere
+
+This page shows how to use Decapodes tooling to explore the harmonics of a discrete manifold. This isn't using any Decapodes specific code, but it 
+is emblematic of a more advanced analysis you might want to do on your Decapode.
+
+In this case we are trying to visualize the roots of the Laplacian on a discrete manifold.
+
+Load the dependencies
+
+```@example Harmonics
+# Meshing:
+using CombinatorialSpaces
+using CoordRefSystems
+using GeometryBasics: Point3
+const Point3D = Point3{Float64};
+
+# Visualization:
+using CairoMakie
+
+# Simulation:
+using LinearAlgebra
+```
+
+Load the mesh
+
+```@example Harmonics
+const RADIUS = 1.0
+s = loadmesh(Icosphere(3, RADIUS));
+sd = EmbeddedDeltaDualComplex2D{Bool, Float64, Point3D}(s);
+subdivide_duals!(sd, Barycenter());
+```
+
+Compute the Laplacian eigenvectors using [LinearAlgebra.eigen](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.eigen). This requires making the sparse Laplacian matrix dense with `collect`. Alternatively, use [Arpack.jl](https://arpack.julialinearalgebra.org/stable/).
+
+```@example Harmonics
+Δ0 = -Δ(0,sd)
+λ = eigen(collect(Δ0))
+```
+
+Let's check that our eigenvalues satisfy the right equation.
+The first eigenvector should be the kernel of the laplacian.
+So the following norm should be close to 0.
+
+```@example Harmonics
+q1 = λ.vectors[:,1]
+norm(Δ0 *q1)
+```
+
+The first eigenvector is boring to visualize, because it is constant.
+So we will make some plots of the second eigenvector.
+If you run this on the desktop, you can use GLMakie and get an interactive plot to explore. We will just draw two angles.
+
+```@example Harmonics
+q = λ.vectors[:,2]
+fig = Figure()
+Label(fig[1, 1, Top()], "Default Angle", padding = (0, 0, 5, 0))
+ax = LScene(fig[1,1], scenekw=(lights=[],))
+msh = CairoMakie.mesh!(ax, s, color=q)
+Colorbar(fig[1,2], msh, size=32)
+# Second Angle
+Label(fig[2, 1, Top()], "Bottom Angle", padding = (0, 0, 5, 0))
+ax = LScene(fig[2,1], scenekw=(lights=[],))
+update_cam!(ax.scene, Vec3f(-1/2,-1/2,1.0/2), Vec3f(1,1,1), Vec3f(0, 0, 1))
+msh = CairoMakie.mesh!(ax, s, color=q)
+Colorbar(fig[2,2], msh, size=32)
+fig
+```
+
+```@example Harmonics
+q = λ.vectors[:,12]
+fig = Figure()
+Label(fig[1, 1, Top()], "Default Angle", padding = (0, 0, 5, 0))
+ax = LScene(fig[1,1], scenekw=(lights=[],))
+msh = CairoMakie.mesh!(ax, s, color=q)
+Colorbar(fig[1,2], msh, size=32)
+# Second Angle
+Label(fig[2, 1, Top()], "Bottom Angle", padding = (0, 0, 5, 0))
+ax = LScene(fig[2,1], scenekw=(lights=[],))
+update_cam!(ax.scene, Vec3f(-1/2,-1/2,1.0/2), Vec3f(1,1,1), Vec3f(0, 0, 1))
+msh = CairoMakie.mesh!(ax, s, color=q)
+Colorbar(fig[2,2], msh, size=32)
+fig
+```
+
+# Exploring solutions with Krylov methods
+
+We can also use the information about the eigenvectors for spectral techniques in solving the equations. Krylov methods are a bridge between
+linear solvers and spectral information.
+
+```julia
+using Krylov
+b = zeros(nv(sd))
+b[1] = 1
+b[end] = -1
+x, stats = Krylov.gmres(Δ0, b, randn(nv(sd)), restart=true, memory=20, atol = 1e-10, rtol=1e-8, history=true, itmax=10000)
+x̂ = x .- sum(x)./length(x)
+norm(x̂)
+stats
+norm(Δ0*(x) - b)
+```
diff --git a/docs/src/index.md b/docs/src/index.md
index 9275ccf0..11a07b1b 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -7,7 +7,6 @@ Decapodes.jl is the synthesis of Applied Category Theory (ACT) techniques for fo
 
 By combining the power of ACT and the DEC, we seek to improve the scientific computing workflow. Decapodes simulations are [hierarchically composable](https://algebraicjulia.github.io/Decapodes.jl/dev/budyko_sellers_halfar/), generalize over [any type of manifold](https://algebraicjulia.github.io/Decapodes.jl/dev/ice_dynamics/), and are [performant and accurate](https://www.cise.ufl.edu/~luke.morris/cism.html) with a declarative domain specific language (DSL) that is [human-readable](https://algebraicjulia.github.io/Decapodes.jl/dev/klausmeier/#Model-Representation).
 
-![Grigoriev Ice Cap Dynamics](grigoriev/grigoriev.gif)
 
 ## Getting started