diff --git a/docs/figs/ch1-07.jl b/docs/figs/ch1-07.jl
index 2e6d5586..ace10843 100644
--- a/docs/figs/ch1-07.jl
+++ b/docs/figs/ch1-07.jl
@@ -12,6 +12,9 @@ using SimpleUnPack
 using Plots
 Plots.default(linewidth=2)
 
+# PNG output in Literate.jl
+PNG(fig) = display("image/png", fig)
+
 # Convenience functions
 hil(x, k) = x / (x + k)
 hil(x, k, n) = hil(x^n, k^n)
@@ -49,3 +52,5 @@ prob = ODEProblem(collins!, u0, tend, ps)
 sol = solve(prob, callback=events)
 
 fig = plot(sol, legend=:right, xlabel = "Time", ylabel="Concentration", title="Fig 1.7")
+
+fig |> PNG
diff --git a/docs/figs/ch1-09.jl b/docs/figs/ch1-09.jl
index 114fe64a..4f683e09 100644
--- a/docs/figs/ch1-09.jl
+++ b/docs/figs/ch1-09.jl
@@ -11,6 +11,9 @@ using SimpleUnPack
 using Plots
 Plots.default(linewidth=2)
 
+# PNG output in Literate.jl
+PNG(fig) = display("image/png", fig)
+
 # Convenience functions
 hil(x, k) = x / (x + k)
 hil(x, k, n) = hil(x^n, k^n)
@@ -64,3 +67,5 @@ p1 = plot(sol, idxs=[:v], ylabel="Membrane potential (mV)", xlabel="", legend=fa
 p2 = plot(sol, idxs = [:m, :h, :n], xlabel="")
 p3 = plot(_istim, sol.t, xlabel="Time (ms)", ylabel="Current", labels="Stimulation current")
 fig = plot(p1, p2, p3, layout=(3, 1), size=(600, 900), leftmargin=5*Plots.mm)
+
+fig |> PNG
diff --git a/docs/figs/ch2-04.jl b/docs/figs/ch2-04.jl
index 4dabdb0a..259fb1b5 100644
--- a/docs/figs/ch2-04.jl
+++ b/docs/figs/ch2-04.jl
@@ -6,9 +6,13 @@ Exponential decay
 using Plots
 Plots.default(linewidth=2)
 
+# PNG output in Literate.jl
+PNG(fig) = display("image/png", fig)
+
+#---
 fig = plot(title= "Fig 2.4")
 for k in 1:3
     plot!(fig, t -> 3 * exp(-k*t), 0., 5., label = "exp(-$(k)t)")
 end
 
-plot!(fig, xlim = (0, 5), ylim=(0, 3.2), xlabel="Time", ylabel="Concentration")
+plot!(fig, xlim = (0, 5), ylim=(0, 3.2), xlabel="Time", ylabel="Concentration") |> PNG
diff --git a/docs/figs/ch2-09.jl b/docs/figs/ch2-09.jl
index 5a3b6faf..abd2923d 100644
--- a/docs/figs/ch2-09.jl
+++ b/docs/figs/ch2-09.jl
@@ -37,3 +37,5 @@ fig = plot(sol, legend=:bottomright, title="Fig 2.9",
     xlims=(0., 4.), ylims=(0., 1.),
     xlabel="Time (sec)", ylabel="Concentration (mM)"
 )
+
+fig |> PNG
diff --git a/docs/figs/ch2-11.jl b/docs/figs/ch2-11.jl
index b6d3d834..8f41c8bb 100644
--- a/docs/figs/ch2-11.jl
+++ b/docs/figs/ch2-11.jl
@@ -10,6 +10,11 @@ using ModelingToolkit
 using Plots
 Plots.default(linewidth=2)
 
+# PNG output in Literate.jl
+PNG(fig) = display("image/png", fig)
+
+#---
+
 rn211 = @reaction_network begin
     k0, 0 --> A
     (k1, km1), A <--> B
@@ -29,7 +34,7 @@ fig = plot(
     xlabel="Time (AU)",
     ylabel="Concentration (AU)",
     title="Fig. 2.11 (Full model)"
-)
+) |> PNG
 
 
 # ## Figure 2.12 : Rapid equilibrium assumption
@@ -80,6 +85,8 @@ plot!(fig,
     ylabel="Concentration (AU)"
 )
 
+fig |> PNG
+
 #===
 ## Figure 2.13: Rapid equilibrium (take 2)
 
@@ -101,6 +108,8 @@ plot!(fig,
     ylabel="Concentration (AU)"
 )
 
+fig |> PNG
+
 
 #===
 ## Figure 2.14 : QSSA
@@ -136,3 +145,5 @@ plot!(fig,
     title="Figure 2.14: Ref vs QSSA",
     xlims=(0.0, tend)
 )
+
+fig |> PNG
diff --git a/docs/figs/ch2-p246.jl b/docs/figs/ch2-p246.jl
index 959b2584..b1b47df8 100644
--- a/docs/figs/ch2-p246.jl
+++ b/docs/figs/ch2-p246.jl
@@ -4,4 +4,10 @@ using DifferentialEquations
 using Plots
 Plots.default(linewidth=2)
 
+# PNG output in Literate.jl
+PNG(fig) = display("image/png", fig)
+
+#---
 fig = ODEProblem((u, p, t) -> p * (1. - u), 0., 10., 1.) |> solve |> plot
+
+fig |> PNG
diff --git a/docs/figs/ch3-03.jl b/docs/figs/ch3-03.jl
index 426d6c7f..67f8796b 100644
--- a/docs/figs/ch3-03.jl
+++ b/docs/figs/ch3-03.jl
@@ -9,6 +9,9 @@ using Catalyst
 using Plots
 Plots.default(linewidth=2)
 
+# PNG output in Literate.jl
+PNG(fig) = display("image/png", fig)
+
 # Reaction neetwork
 rn = @reaction_network begin
     (k1, km1), S + E <--> ES
@@ -42,6 +45,8 @@ sol = solve(prob)
 @unpack S, ES, E, P = osys
 fig = plot(sol, idxs=[S, ES, E, P], xlabel="Time (AU)", ylabel="Concentration (AU)", legend=:right, title="Fig 3.03")
 
+fig |> PNG
+
 #---
 rn303mm = @reaction_network begin
     mm(S, k2 * ET, (km1 + k2) / k1), S ⇒ P ## using \Rightarrow
@@ -72,3 +77,5 @@ plot!(fig, title="Fig. 3.03",
     xlabel="Time (AU)", ylabel="Concentration (AU)",
     xlims=(0., tend), ylims=(0., 5.), legend=:right
 )
+
+fig |> PNG
diff --git a/docs/figs/ch3-13.jl b/docs/figs/ch3-13.jl
index 98978736..8f36fe6c 100644
--- a/docs/figs/ch3-13.jl
+++ b/docs/figs/ch3-13.jl
@@ -7,6 +7,10 @@ Generalized mass action (GMA) vs. Michaelis-Menten rate laws
 using Plots
 Plots.default(linewidth=2)
 
+# PNG output in Literate.jl
+PNG(fig) = display("image/png", fig)
+
+#---
 mm = t -> 2t / (1+t)
 gma = t -> t^0.4
 
@@ -16,3 +20,5 @@ fig = plot(
     xlabel= "Substrate concentration (AU)",
     ylabel="Reaction rate (AU)"
 )
+
+fig |> PNG
diff --git a/docs/figs/ch4-01.jl b/docs/figs/ch4-01.jl
index 97d992f0..983758e5 100644
--- a/docs/figs/ch4-01.jl
+++ b/docs/figs/ch4-01.jl
@@ -10,8 +10,10 @@ using SimpleUnPack
 using Plots
 Plots.default(linewidth=2)
 
-# Convenience functions
+# PNG output in Literate.jl
+PNG(fig) = display("image/png", fig)
 
+# Convenience functions
 hil(x, k) = x / (x + k)
 hil(x, k, n) = hil(x^n, k^n)
 
@@ -54,6 +56,8 @@ sols = [solve(ODEProblem(model41!, u0, tend, ps1)) for u0 in u0s];
 #---
 fig = plot(sols[1], xlabel="Time", ylabel="Concentration", title="Fig. 4.2 A (Time series)", labels=["[A]" "[B]"])
 
+fig |> PNG
+
 # ## Fig. 4.2 B (Phase plot)
 
 fig = plot( sols[1], idxs=(:a, :b),
@@ -63,6 +67,8 @@ fig = plot( sols[1], idxs=(:a, :b),
     ylims=(0.0, 2.0), xlims=(0.0, 2.0)
 )
 
+fig |> PNG
+
 # ## Fig. 4.3 A (Multiple time series)
 
 fig = plot(title="Fig. 4.3A (Multiple time series)")
@@ -73,6 +79,8 @@ end
 
 plot!(fig, xlabel="Time", ylabel="Concentration")
 
+fig |> PNG
+
 # ## Fig. 4.3 B (Phase plot)
 
 fig = plot(title="Fig. 4.3 B (Phase plot)")
@@ -83,6 +91,8 @@ end
 
 plot!(fig, xlabel="[A]", ylabel="[B]", xlims=(0., 2.), ylims=(0., 2.), aspect_ratio=:equal)
 
+fig |> PNG
+
 # Let's sketch vector fields in phase plots.
 ∂A = function (x, y)
     ∂A41((x, y), ps1, 0.0)
@@ -115,6 +125,8 @@ end
 
 plot!(fig, size=(600, 600), xlims=(rxy[1], rxy[end]), ylims=(rxy[1], rxy[end]), xlabel="[A]", ylabel="[B]")
 
+fig |> PNG
+
 # ## Figure 4.5A
 
 fig = plot(title="Fig. 4.5 A (Phase plot with nullclines)")
@@ -131,6 +143,8 @@ contour!(fig, 0:0.01:2, 0:0.01:2, ∂B, levels=[0], cbar=false, line=(:black, :d
 plot!(fig, Float64[], Float64[], line=(:black, :dash), label="B nullcline") ## Adding a fake line for the legend of B nullcline
 plot!(fig, xlim=(0., 2.), ylim=(0., 2.), legend=:bottomright, size=(600, 600), xlabel="[A]", ylabel="[B]", aspect_ratio=:equal)
 
+fig |> PNG
+
 # ## Figure 4.5 B
 
 ## Vector field
@@ -143,3 +157,5 @@ plot!(fig, Float64[], Float64[], line=(:black), label="A nullcline")  ## Adding
 contour!(fig, 0:0.01:2, 0:0.01:2, ∂B, levels=[0], cbar=false, line=(:black, :dash))
 plot!(fig, Float64[], Float64[], line=(:black, :dash), label="B nullcline") ## Adding a fake line for the legend of B nullcline
 plot!(fig, xlim=(0., 2.), ylim=(0., 2.), legend=:bottomright, size=(600, 600), xlabel="[A]", ylabel="[B]")
+
+fig |> PNG
diff --git a/docs/figs/ch4-07.jl b/docs/figs/ch4-07.jl
index 401aeb3e..fa9c7441 100644
--- a/docs/figs/ch4-07.jl
+++ b/docs/figs/ch4-07.jl
@@ -7,6 +7,9 @@ using SimpleUnPack
 using Plots
 Plots.default(linewidth=2)
 
+# PNG output in Literate.jl
+PNG(fig) = display("image/png", fig)
+
 # Convenience functions
 hil(x, k) = x / (x + k)
 hil(x, k, n) = hil(x^n, k^n)
@@ -46,6 +49,8 @@ ax1 = plot(sol1, xlabel="Time", ylabel="Concentration", legend=:right, title= "F
 ax2 = plot(sol2, xlabel="Time", ylabel="Concentration", legend=:right, title= "Fig 4.7A (2)")
 fig = plot(ax1, ax2, layout=(2, 1), size=(600, 600))
 
+fig |> PNG
+
 # ## Fig 4.7 B
 
 ∂F47 = function (x, y; scale=20)
@@ -77,6 +82,8 @@ contour!(fig, 0:0.01:5, 0:0.01:5, ∂B, levels=[0], cbar=false, line=(:black, :d
 plot!(fig, Float64[], Float64[], line=(:black, :dash), label="B nullcline")
 plot!(fig, xlim=(0, 5), ylim=(0, 5), legend=:topright, size=(600, 600), xlabel="[A]", ylabel="[B]")
 
+fig |> PNG
+
 #===
 ## Fig 4.8
 
@@ -93,6 +100,8 @@ ax1 = plot(sol1, xlabel="Time", ylabel="Concentration", legend=:right, title= "F
 ax2 = plot(sol2, xlabel="Time", ylabel="Concentration", legend=:right, title= "Fig 4.8A (2)")
 fig = plot(ax1, ax2, layout=(2, 1), size=(600, 600))
 
+fig |> PNG
+
 #---
 ∂F48 = function (x, y; scale=20)
     da = ∂A47((x, y), ps2, 0.0)
@@ -119,6 +128,8 @@ contour!(fig, 0:0.01:5, 0:0.01:5, ∂B, levels=[0], cbar=false, line=(:black, :d
 plot!(fig, Float64[], Float64[], line=(:black, :dash), label="B nullcline")
 plot!(fig, xlim=(0, 5), ylim=(0, 5), legend=:topright, size=(600, 600), xlabel="[A]", ylabel="[B]")
 
+fig |> PNG
+
 #===
 ## Fig 4.8 C
 around the unstable steady-state
@@ -137,6 +148,8 @@ contour!(fig, 1:0.01:1.5, 1:0.01:1.5, ∂B, levels=[0], cbar=false, line=(:black
 plot!(fig, Float64[], Float64[], line=(:black, :dash), label="B nullcline")
 plot!(fig, xlim=(1, 1.5), ylim=(1, 1.5), legend=:topright, size=(600, 600), xlabel="[A]", ylabel="[B]")
 
+fig |> PNG
+
 # Another way to draw nullclines is to find the analytical solution when dA (or dB) is zero. And then sketch the nullclines in a parameteric plot.
 
 nca47(b, p) = p.k1 / p.k3 * hil(1, b, p.n1)
@@ -151,3 +164,5 @@ pls = map((8.0, 16.0, 20.0, 35.0)) do k1
 end
 
 fig = plot(pls..., size = (800, 800))
+
+fig |> PNG
diff --git a/docs/figs/ch4-11.jl b/docs/figs/ch4-11.jl
index 88c419e1..fab3f380 100644
--- a/docs/figs/ch4-11.jl
+++ b/docs/figs/ch4-11.jl
@@ -7,6 +7,10 @@ Surface plots reference: [surface plots @ PlotsGallery.jl](https://goropikari.gi
 using Plots
 Plots.default(linewidth=2)
 
+# PNG output in Literate.jl
+PNG(fig) = display("image/png", fig)
+
+#---
 z1(x, y) = x^2 + 0.5y^2
 z2(x, y) = (.2x^2-1)^2 + y^2
 x1 = y1 = range(-1.0, 1.0, 41)
@@ -18,3 +22,5 @@ p3 = surface(x2, y2, z2, title="Double-well potential")
 p4 = contourf(x2, y2, z2)
 
 fig = plot(p1, p2, p3, p4, size=(800, 600))
+
+fig |> PNG
diff --git a/docs/figs/ch4-15.jl b/docs/figs/ch4-15.jl
index 6a8c0bc1..a1ab4ae0 100644
--- a/docs/figs/ch4-15.jl
+++ b/docs/figs/ch4-15.jl
@@ -12,6 +12,9 @@ using SimpleUnPack
 using Plots
 Plots.default(linewidth=2)
 
+# PNG output in Literate.jl
+PNG(fig) = display("image/png", fig)
+
 #---
 
 function dA415(u, p, t)
@@ -48,6 +51,8 @@ end
 
 fig = plot(sols[1], xlabel="Time", ylabel="Concentration", title ="Fig 4.15 (A)", xlims=(0., 8.))
 
+fig |> PNG
+
 # ## Fig 4.15 (B)
 
 ∂F415 = function (x, y; scale=20)
@@ -82,6 +87,8 @@ for sol in sols
 end
 plot!(fig, xlim=(0, 4), ylim=(0, 4), legend=:topright, size=(600, 600), xlabel="[A]", ylabel="[B]")
 
+fig |> PNG
+
 #===
 ## Fig 4.16 A
 
@@ -103,6 +110,8 @@ end;
 
 fig = plot(sols[1], xlabel="Time", ylabel="Concentration", title ="Fig 4.16(A)", xlims=(0., 8.))
 
+fig |> PNG
+
 # ## Fig 4.16 b
 ∂F416 = function (x, y; scale=20)
     da = dA415((x, y), ps2, 0.0)
@@ -135,6 +144,8 @@ for sol in sols
 end
 plot!(fig, xlim=(0, 4), ylim=(0, 4), legend=:topright, size=(600, 600), xlabel="[A]", ylabel="[B]")
 
+fig |> PNG
+
 # ## Fig 4.17
 
 sol = solve(ODEProblem(model415!, LVector(a=2.0, b=1.5), 10.0, ps2))
@@ -147,3 +158,5 @@ plot!(fig, identity, 0, 0, line=(:black), label="A nullcline")
 contour!(fig, 1:0.01:3, 1:0.01:3, ∂B, levels=[0], cbar=false, line=(:black, :dash))
 plot!(fig, identity, 0, 0, line=(:black, :dash), label="B nullcline")
 plot!(fig, xlims=(1, 3), ylims=(1, 3), legend=:topright, size=(600, 600), xlabel="[A]", ylabel="[B]")
+
+fig |> PNG
diff --git a/docs/figs/ch4-18.jl b/docs/figs/ch4-18.jl
index 82e20c89..b9ecdbde 100644
--- a/docs/figs/ch4-18.jl
+++ b/docs/figs/ch4-18.jl
@@ -12,6 +12,9 @@ using LabelledArrays
 using Plots
 Plots.default(linewidth=2)
 
+# PNG output in Literate.jl
+PNG(fig) = display("image/png", fig)
+
 #---
 hil(x, k) = x / (x + k)
 hil(x, k, n) = hil(x^n, k^n)
@@ -39,3 +42,5 @@ fig = plot(
     xlabel = "K1" , ylabel= "Steady state [A]",
     legend=nothing, ylim=(0, 4), xlim=(0, 1000)
 )
+
+fig |> PNG
diff --git a/docs/figs/ch4-22.jl b/docs/figs/ch4-22.jl
index 5cdcb292..e59afad6 100644
--- a/docs/figs/ch4-22.jl
+++ b/docs/figs/ch4-22.jl
@@ -10,6 +10,9 @@ using ForwardDiff
 using Plots
 Plots.default(linewidth=2)
 
+# PNG output in Literate.jl
+PNG(fig) = display("image/png", fig)
+
 # The curve
 f = x -> 3 / (x-2)
 
@@ -26,3 +29,5 @@ plot!(fig, g, 2.7, 5.3, lab="Tangent line")
 plot!(fig, xlabel="Reaction rate", ylabel="Inhibitor concentration",
       xlims=(2., 8.), ylims=(0., 4.)
 )
+
+fig |> PNG
diff --git a/docs/figs/ch5-10.jl b/docs/figs/ch5-10.jl
index e0bca881..e0a8dce6 100644
--- a/docs/figs/ch5-10.jl
+++ b/docs/figs/ch5-10.jl
@@ -10,6 +10,9 @@ using SimpleUnPack
 using Plots
 Plots.default(linewidth=2)
 
+# PNG output in Literate.jl
+PNG(fig) = display("image/png", fig)
+
 #---
 hil(x, k=one(x)) = x / (x + k)
 hil(x, k, n) = hil(x^n, k^n)
@@ -62,6 +65,8 @@ sol = solve(prob)
 
 fig = plot(sol, title="Figure 5.10", xlabel="Time (hr)", ylabel="Concentration (μM)", xlims=(0, 1), legend=:right)
 
+fig |> PNG
+
 # ## Figure 5.11 A
 
 rx = range(0, 1200, 101)
@@ -101,6 +106,8 @@ end
 
 plot!(fig, xlims=(0, 1200), ylims=(0, 6), xlabel="AdoMet (μM)", ylabel="AdoHcy (μM)", legend=:bottomright)
 
+fig |> PNG
+
 # ## Figure 5.11 B
 # Increase methionine level
 
@@ -139,3 +146,5 @@ for sol in sols
 end
 
 plot!(fig, xlims=(0, 1200), ylims=(0, 6), xlabel="AdoMet (μM)", ylabel="AdoHcy (μM)", legend=:bottomright)
+
+fig |> PNG
diff --git a/docs/figs/ch6-03.jl b/docs/figs/ch6-03.jl
index e1d9c07e..faaa2861 100644
--- a/docs/figs/ch6-03.jl
+++ b/docs/figs/ch6-03.jl
@@ -9,6 +9,9 @@ using DifferentialEquations
 using Plots
 Plots.default(linewidth=2)
 
+# PNG output in Literate.jl
+PNG(fig) = display("image/png", fig)
+
 #---
 rn = @reaction_network begin
     (k1 * L, km1), R <--> RL
@@ -52,6 +55,8 @@ fig = plot(sol, idxs=[RL, Ps], labels= ["RL" "P*"])
 plot!(fig, t -> 3 * (1<=t<=3), label="Ligand", line=(:black, :dash), linealpha=0.7)
 plot!(fig, title="Fig. 6.3 (A)", xlabel="Time", ylabel="Concentration")
 
+fig |> PNG
+
 # ## Fig 6.3 B
 
 lrange = 0:0.01:1
@@ -70,3 +75,5 @@ pstar = map(s->s[Ps], sim)
 rl = map(s->s[RL], sim)
 fig = plot(lrange, [pstar rl], label=["P*" "RL"], title="Fig. 6.3 (B)",
 xlabel="Ligand", ylabel="Steady-state concentration", xlims=(0, 1), ylims=(0, 8))
+
+fig |> PNG
diff --git a/docs/figs/ch6-05.jl b/docs/figs/ch6-05.jl
index 79cebe1e..ef3a1599 100644
--- a/docs/figs/ch6-05.jl
+++ b/docs/figs/ch6-05.jl
@@ -10,8 +10,10 @@ using DifferentialEquations
 using Plots
 Plots.default(linewidth=2)
 
-#---
+# PNG output in Literate.jl
+PNG(fig) = display("image/png", fig)
 
+#---
 rn = @reaction_network begin
     (kRL * L, kRLm), R <--> RL
     kGa, G + RL --> Ga + Gbg + RL
@@ -56,6 +58,8 @@ sol = solve(prob, Rodas5(), callback=cbs, abstol=1e-8, reltol=1e-8, saveat=0:10:
 @unpack RL, Ga = osys
 fig = plot(sol, idxs=[RL, Ga], title="Fig 6.05 (A)", xlabel="Time", ylabel="Abundance")
 
+fig |> PNG
+
 # ## Fig 6.5 B
 
 lrange = range(0, 20 * 1e-9, 101)
@@ -74,3 +78,5 @@ ga = map(s->s[Ga], sim)
 rl = map(s->s[RL], sim)
 fig = plot(lrange .* 1e9, [ga rl], label=["Ga" "RL"], title="Fig. 6.5 (B)",
 xlabel="Ligand (nM)", ylabel="Steady-state abundance", xlims=(0, 20), ylims=(0, 3500))
+
+fig |> PNG
diff --git a/docs/figs/ch6-07.jl b/docs/figs/ch6-07.jl
index 8299c891..1746bc3b 100644
--- a/docs/figs/ch6-07.jl
+++ b/docs/figs/ch6-07.jl
@@ -7,9 +7,14 @@ Goldbeter Koshland switch
 using Plots
 Plots.default(linewidth=2)
 
+# PNG output in Literate.jl
+PNG(fig) = display("image/png", fig)
+
 f607(w, K1, K2) = w * (1 - w + K1)/((1 - w) * (w + K2))
 
 fig = plot(title = "Fig 6.07", xlims=(0, 3), ylims=(0, 1))
 plot!(fig, x -> f607(x, 1, 1), y -> y, 0, 1, label="K1=K2=1")
 plot!(fig, x -> f607(x, 0.1, 0.1), y -> y, 0, 1, label="K1=K2=0.1")
 plot!(fig, x -> f607(x, 0.01, 0.01), y -> y, 0, 1, label="K1=K2=0.01")
+
+fig |> PNG
diff --git a/docs/figs/ch6-14.jl b/docs/figs/ch6-14.jl
index 1b2ddb11..15e8b92e 100644
--- a/docs/figs/ch6-14.jl
+++ b/docs/figs/ch6-14.jl
@@ -10,6 +10,9 @@ using DifferentialEquations
 using Plots
 Plots.default(linewidth=2)
 
+# PNG output in Literate.jl
+PNG(fig) = display("image/png", fig)
+
 #---
 
 rn = @reaction_network begin
@@ -71,3 +74,5 @@ sol = solve(prob, callback=cbs)
 
 @unpack Am = osys
 fig = plot(sol, idxs=[Am], title="Fig 6.14", xlabel="Time", ylabel="Active CheA ([Am])", ylims=(0.01, 0.04), legend=false)
+
+fig |> PNG
diff --git a/docs/figs/ch6-16.jl b/docs/figs/ch6-16.jl
index ad8b3eb0..405f70e5 100644
--- a/docs/figs/ch6-16.jl
+++ b/docs/figs/ch6-16.jl
@@ -10,6 +10,9 @@ using DifferentialEquations
 using Plots
 Plots.default(linewidth=2)
 
+# PNG output in Literate.jl
+PNG(fig) = display("image/png", fig)
+
 #---
 rn = @reaction_network begin
     (k1, k2), 0 <--> C8
@@ -77,3 +80,5 @@ sol = solve(prob, callback=cbs)
 
 @unpack C8s, C3s = osys
 fig = plot(sol, idxs=[C8s, C3s], title="Fig 6.16", xlabel="Time", ylabel="Concentration", legend=:right, rightmargin=5*Plots.mm)
+
+fig |> PNG
diff --git a/docs/figs/ch6-18.jl b/docs/figs/ch6-18.jl
index 9c9e7083..227bca96 100644
--- a/docs/figs/ch6-18.jl
+++ b/docs/figs/ch6-18.jl
@@ -10,6 +10,9 @@ using DifferentialEquations
 using Plots
 Plots.default(linewidth=2)
 
+# PNG output in Literate.jl
+PNG(fig) = display("image/png", fig)
+
 #---
 rn = @reaction_network begin
     (k1 * I, km1), R <--> RI
@@ -52,6 +55,8 @@ sol = solve(prob)
 @unpack C, RIC, RICC = osys
 fig = plot(sol, idxs=[C, RIC, RICC], title="Fig 6.18 (A)", xlabel="Time", ylabel="Abundance", legend=:topright)
 
+fig |> PNG
+
 # ## Fig 6.18 (B)
 
 idx = findfirst(isequal(I), parameters(osys))
@@ -65,3 +70,5 @@ sol = solve(prob, callback=cbs)
 
 @unpack C = osys
 fig = plot(sol, idxs=[C], title="Fig 6.18 (B)", xlabel="Time", ylabel="Ca concentration", legend=false, ylim=(0, 2.5))
+
+fig |> PNG
diff --git a/docs/figs/ch7-07.jl b/docs/figs/ch7-07.jl
index 787c4133..bdb5a862 100644
--- a/docs/figs/ch7-07.jl
+++ b/docs/figs/ch7-07.jl
@@ -10,6 +10,9 @@ using DifferentialEquations
 using Plots
 Plots.default(linewidth=2)
 
+# PNG output in Literate.jl
+PNG(fig) = display("image/png", fig)
+
 #---
 rn = @reaction_network begin
     a1 / (1+RToverK1*(K2/(K2+L))^4), 0 --> M
@@ -75,6 +78,8 @@ end
 
 plot!(fig, lac, 0, 2500, label="External lactose (μM)")
 
+fig |> PNG
+
 # ## Fig 7.07 (B)
 # Compare the original model and the modified model
 
@@ -141,3 +146,5 @@ fig = plot(lerange, [sim sim_mod], label=["Original" "Modified"],
     ylabel="β-galactosidase",
     title="Fig 7.7 (B)"
 )
+
+fig |> PNG
diff --git a/docs/figs/ch7-11.jl b/docs/figs/ch7-11.jl
index cc8b24c7..02f352df 100644
--- a/docs/figs/ch7-11.jl
+++ b/docs/figs/ch7-11.jl
@@ -10,6 +10,9 @@ using SimpleUnPack
 using Plots
 Plots.default(linewidth=2)
 
+# PNG output in Literate.jl
+PNG(fig) = display("image/png", fig)
+
 #---
 
 function model711(u, p, t)
@@ -61,6 +64,8 @@ quiver!(fig, xx, yy, quiver=∂F, line=(:lightblue), arrow=(:closed))
 
 plot!(fig, xlims=(0, 250), ylims=(0, 250), xlabel="[cI] (nM)", ylabel="[cro] (nM)", aspect_ratio=:equal, legend=:top, size=(600, 600))
 
+fig |> PNG
+
 # ## Fig 7.11 (B)
 
 ps2 = merge(ps1, (; delta_r=ps1.delta_r * 10))
@@ -88,3 +93,5 @@ plot!(fig, Float64[], Float64[], line=(:black, :dash), label="C nullcline")
 quiver!(fig, xx, yy, quiver=∂F, line=(:lightblue), arrow=(:closed))
 
 plot!(fig, xlims=(0, 250), ylims=(0, 250), xlabel="[cI] (nM)", ylabel="[cro] (nM)", aspect_ratio=:equal, legend=:top, size=(600, 600))
+
+fig |> PNG
diff --git a/docs/figs/ch7-17.jl b/docs/figs/ch7-17.jl
index 3c595e54..d4eb65b1 100644
--- a/docs/figs/ch7-17.jl
+++ b/docs/figs/ch7-17.jl
@@ -10,6 +10,9 @@ using DifferentialEquations
 using Plots
 Plots.default(linewidth=2)
 
+# PNG output in Literate.jl
+PNG(fig) = display("image/png", fig)
+
 #---
 rn = @reaction_network begin
     (a / (k^n + Z^n), b), 0 <--> X
@@ -41,5 +44,9 @@ sol = solve(prob)
 #---
 fig = plot(sol, title="Fig 7.17 (A)", xlabel="Time", ylabel="Concentration")
 
+fig |> PNG
+
 #---
 fig = plot(sol, idxs=(1, 2, 3), title="Fig 7.17 (B)", legend=false, size=(600, 600))
+
+fig |> PNG
diff --git a/docs/hw-01.jl b/docs/hw-01.jl
index 1a880e04..fbdf3b75 100644
--- a/docs/hw-01.jl
+++ b/docs/hw-01.jl
@@ -66,8 +66,13 @@ analytical(t) = 1 - exp(-t)
 using Plots
 Plots.default(linewidth=2)
 
-plot(sol.t, sol.u, label="FE method")
-plot!(analytical, sol.t[begin], sol.t[end], label = "Analytical solution", linestyle=:dash)
+# PNG output in Literate.jl
+PNG(fig) = display("image/png", fig)
+
+fig = plot(sol.t, sol.u, label="FE method")
+plot!(fig, analytical, sol.t[begin], sol.t[end], label = "Analytical solution", linestyle=:dash)
+
+fig |> PNG
 
 #===
 
diff --git a/docs/intro-02-plotting.jl b/docs/intro-02-plotting.jl
index a5227e07..23816230 100644
--- a/docs/intro-02-plotting.jl
+++ b/docs/intro-02-plotting.jl
@@ -9,6 +9,8 @@
 ===#
 
 using Plots
+# PNG output in Literate.jl
+PNG(fig) = display("image/png", fig)
 
 #---
 f(x) = sin(sin(x) + 1)
@@ -114,5 +116,6 @@ quiver(xx, yy, quiver=∇f, aspect_ratio=:equal, line=(:black), arrow=(:closed))
 
 # Save the current figure
 savefig("vfield.png")
+
 # Save the figure saved in a variable
 savefig(fig, "subplots.png")
diff --git a/docs/intro-03-diffeq.jl b/docs/intro-03-diffeq.jl
index ab859374..2f44d772 100644
--- a/docs/intro-03-diffeq.jl
+++ b/docs/intro-03-diffeq.jl
@@ -22,6 +22,9 @@ using DifferentialEquations
 using Plots
 Plots.default(linewidth=2)
 
+# PNG output in Literate.jl
+PNG(fig) = display("image/png", fig)
+
 #===
 
 ### Exponential decay model
@@ -55,7 +58,9 @@ prob = ODEProblem(expdecay, u0, tspan, p)
 sol = solve(prob)
 
 # Visualize the solution
-plot(sol, legend=:right)
+fig = plot(sol, legend=:right)
+
+fig |> PNG
 
 # Solution at time t=1.0 (with interpolation)
 sol(1.0)
@@ -117,7 +122,9 @@ prob = ODEProblem(sir!, u0, tspan, p)
 sol = solve(prob)
 
 # Visualize the solution
-plot(sol, label=["S" "I" "R"], legend=:right)
+fig = plot(sol, label=["S" "I" "R"], legend=:right)
+
+fig |> PNG
 
 #===
 ### Lorenz system
@@ -155,15 +162,19 @@ prob = ODEProblem(lorenz!,u0,tspan,p)
 sol = solve(prob)
 
 # x-y-z time-series
-plot(sol)
+fig = plot(sol)
+
+fig |> PNG
 
 # `idxs=(1, 2, 3)` makes a phase plot with 1st, 2nd, and the 3rd state variable. With `LabelledArrays`, you can use symbols instead of index numbers.
 
-plot(sol, idxs=(:x, :y, :z))
+fig = plot(sol, idxs=(:x, :y, :z))
+fig |> PNG
 
 # The zeroth variable in `idxs` is the independent variable (usually time). The below command plots the time series of the second state variable (`y`).
 
-plot(sol, idxs=(0, 2))
+fig = plot(sol, idxs=(0, 2))
+fig |> PNG
 
 # ## Saving simulation results
 
@@ -206,7 +217,8 @@ tspan = (0.0, 2.0)
 prob = ODEProblem(expdecaySys, u0, tspan, p)
 sol = solve(prob)
 
-plot(sol)
+fig = plot(sol)
+fig |> PNG
 
 # ### SIR model
 
@@ -229,7 +241,8 @@ tspan = (0.0, 20.0)
 prob = ODEProblem(sirSys, u0, tspan, p)
 sol = solve(prob)
 
-plot(sol)
+fig = plot(sol)
+fig |> PNG
 
 #===
 ## Using Catalyst.jl for chemical reaction networks
@@ -255,7 +268,8 @@ tspan = (0., 2.)
 prob = ODEProblem(decay_rn, u0, tspan, p)
 sol = solve(prob)
 
-plot(sol, title="Exponential Decay")
+fig = plot(sol, title="Exponential Decay")
+fig |> PNG
 
 # ### SIR model
 
@@ -271,4 +285,5 @@ tspan = (0., 20.)
 prob = ODEProblem(sir_rn, u0, tspan, p)
 sol = solve(prob)
 
-plot(sol, legend=:right, title = "SIR Model")
+fig = plot(sol, legend=:right, title = "SIR Model")
+fig |> PNG
diff --git a/docs/intro-04-gillespie.jl b/docs/intro-04-gillespie.jl
index 453e345a..f004ae41 100644
--- a/docs/intro-04-gillespie.jl
+++ b/docs/intro-04-gillespie.jl
@@ -12,6 +12,8 @@ using Random            ## randexp()
 Random.seed!(2022)
 Plots.default(fmt=:png)
 
+# PNG output in Literate.jl
+PNG(fig) = display("image/png", fig)
 
 #===
 Stochastic chemical reaction: Gillespie Algorithm (direct method)
@@ -85,15 +87,19 @@ soldirect = ssa_direct(model, u0, tend, parameters, stoich)
 solfirst = ssa_first(model, u0, tend, parameters, stoich)
 
 # Plot the solution from the direct method
-plot(soldirect.t, soldirect.u,
+fig = plot(soldirect.t, soldirect.u,
     xlabel="time", ylabel="# of molecules",
     title = "SSA (direct method)", label=["A" "B"])
 
+fig |> PNG
+
 # Plot the solution by the first reaction method
-plot(solfirst.t, solfirst.u,
+fig = plot(solfirst.t, solfirst.u,
     xlabel="time", ylabel="# of molecules",
     title = "SSA (1st reaction method)", label=["A" "B"])
 
+fig |> PNG
+
 # Running 50 simulations
 numRuns = 50
 
@@ -121,13 +127,13 @@ for sol in sols
     plot!(fig1, sol.t, sol.u, linecolor=[:blue :red], linealpha=0.05, label=false)
 end
 
-fig1
+fig1 |> PNG
 
 # Plot averages
 plot!(fig1, a_avg, 0.0, tend, linecolor=:black, linewidth=3, linestyle = :solid, label="Avarage [A]")
 plot!(fig1, b_avg, 0.0, tend, linecolor=:black, linewidth=3, linestyle = :dash, label="Avarage [B]")
 
-fig1
+fig1 |> PNG
 
 #===
 ## Using Catalyst
@@ -158,19 +164,19 @@ jumpProb = JumpProblem(rn, dprob, Direct())
 sol = solve(jumpProb, SSAStepper())
 fig = plot(sol)
 
-fig
+fig |> PNG
 
 # Parallel ensemble simulation
 ensprob = EnsembleProblem(jumpProb)
 sim = solve(ensprob, SSAStepper(), EnsembleThreads(); trajectories=50)
 fig = plot(sim, alpha=0.1, color=[:blue :red])
 
-fig
+fig |> PNG
 #---
 
 summ = EnsembleSummary(sim, 0:0.1:10)
 fig = plot(summ,fillalpha=0.5)
-fig
+fig |> PNG
 
 #===
 **See also** the [JumpProcesses.jl docs](https://docs.sciml.ai/JumpProcesses/stable/api/#JumpProcesses.ConstantRateJump) about discrete stochastic examples.