diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 6e53c52..9454bb8 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -95,7 +95,7 @@ jobs:
           julia residuals_terms.jl
           julia residuals_terms_px.jl
           julia tau_vs_e.jl
-          #julia waveform_residuals_px.jl
+          julia waveform_residuals_px.jl
           #python simulation_example.py
           #python enterprise_example.py
           cd ..
diff --git a/examples/tau_vs_e.jl b/examples/tau_vs_e.jl
index 97c1673..72b26af 100644
--- a/examples/tau_vs_e.jl
+++ b/examples/tau_vs_e.jl
@@ -10,7 +10,7 @@ es = Eccentricity.(LinRange(0.01, 0.99, 1000))
 τs = τ_from_e.(es)
 
 fig = Figure()
-ax = Axis(fig[1,1])
+ax = Axis(fig[1, 1])
 CairoMakie.lines!([τ.τ for τ in τs], [e.e for e in es])
 ax.xlabel = "τ"
 ax.ylabel = "e"
diff --git a/examples/waveform_residuals_px.jl b/examples/waveform_residuals_px.jl
index aebef20..396d38e 100644
--- a/examples/waveform_residuals_px.jl
+++ b/examples/waveform_residuals_px.jl
@@ -1,7 +1,7 @@
 """Comparison of eccentric waveforms and residuals (+/x polarizations) for different eccentricities."""
 
 using GWecc
-using PyPlot
+using CairoMakie
 
 println("Running ", PROGRAM_FILE)
 
@@ -17,6 +17,8 @@ proj = ProjectionParams(1e-9, 0.0, 1.0, 0.0, 0.0)
 
 ts = Time.(LinRange(0, 10 * year, 5000))
 
+fig = CairoMakie.Figure()
+
 for (idx, e_init) in enumerate(Eccentricity.([0.1, 0.4, 0.8]))
     coeffs = EvolvCoeffs(mass, n_init, e_init)
 
@@ -24,21 +26,21 @@ for (idx, e_init) in enumerate(Eccentricity.([0.1, 0.4, 0.8]))
     hps = [hpx[1] for hpx in hpxs]
     hxs = [hpx[2] for hpx in hpxs]
 
-    subplot(320 + 2 * idx - 1)
-    plot([t.t for t in ts] / year, hps)
-    plot([t.t for t in ts] / year, hxs)
-    ylabel("\$h_{+,\\times}\$")
-    xlabel("t (year)")
+    ax1 = CairoMakie.Axis(fig[idx, 1])
+    CairoMakie.lines!(ax1, [t.t for t in ts] / year, hps)
+    CairoMakie.lines!(ax1, [t.t for t in ts] / year, hxs)
+    ax1.ylabel = "h+, hx"
+    ax1.xlabel = "t (year)"
 
     spxs = [residual_px(mass, coeffs, l0p, proj, false, dt) for dt in ts]
     sps = [spx[1] for spx in spxs]
     sxs = [spx[2] for spx in spxs]
 
-    subplot(320 + 2 * idx)
-    plot([t.t for t in ts] / year, sps)
-    plot([t.t for t in ts] / year, sxs)
-    ylabel("\$s_{+,\\times}\$")
-    xlabel("t (year)")
+    ax2 = CairoMakie.Axis(fig[idx, 2])
+    CairoMakie.lines!(ax2, [t.t for t in ts] / year, sps)
+    CairoMakie.lines!(ax2, [t.t for t in ts] / year, sxs)
+    ax2.ylabel = "s+, sx"
+    ax2.xlabel = "t (year)"
 end
 
-show()
+save("waveform_residuals_px.pdf", fig)