diff --git a/quarto/ODEs/differential_equations.qmd b/quarto/ODEs/differential_equations.qmd
index 5d8d9e4e..992dbb9f 100644
--- a/quarto/ODEs/differential_equations.qmd
+++ b/quarto/ODEs/differential_equations.qmd
@@ -384,15 +384,16 @@ Finally, we note that the `ModelingToolkit` package provides symbolic-numeric co
 
 ```{julia}
 #| hold: true
-@parameters t γ g
-@variables x(t) y(t)
+@parameters γ g
+@variables t x(t) y(t)
 D = Differential(t)
 
 eqs = [D(D(x)) ~ -γ * D(x),
        D(D(y)) ~ -g - γ * D(y)]
 
-@named sys = ODESystem(eqs)
+@named sys = ODESystem(eqs, t, [x,y], [γ,g])
 sys = ode_order_lowering(sys) # turn 2nd order into 1st
+sys = structural_simplify(sys)
 
 u0 = [D(x) => vxy₀[1],
       D(y) => vxy₀[2],
diff --git a/quarto/_common_code.qmd b/quarto/_common_code.qmd
index 45327f90..fcc15752 100644
--- a/quarto/_common_code.qmd
+++ b/quarto/_common_code.qmd
@@ -1,4 +1,5 @@
 ```{julia}
+#| output: false
 #| echo: false
 
 ## Formatting options are included here; not in CalculusWithJulia.WeaveSupport
@@ -7,18 +8,20 @@ nothing
 ```
 
 ```{julia}
+#| output: false
 #| echo: false
 fig_size=(800, 600)
 nothing
 ```
 
 ```{julia}
+#| output: false
 #| echo: false
 
 import Logging
 Logging.disable_logging(Logging.Info) # or e.g. Logging.Info
 Logging.disable_logging(Logging.Warn)
-
+nothing
 ```
 
 ```{julia}
@@ -36,6 +39,7 @@ end
 ```
 
 ```{julia}
+#| output: false
 #| echo: false
 # ImageFile
 ## WeaveSupport from CalculusWithJulia package
@@ -220,4 +224,6 @@ function Base.show(io::IO, m::MIME"text/plain", x::HTMLoutput)
     println(io, caption)
     return nothing
 end
+
+nothing
 ```
diff --git a/quarto/_quarto.yml b/quarto/_quarto.yml
index b16a8f8a..e576fb7b 100644
--- a/quarto/_quarto.yml
+++ b/quarto/_quarto.yml
@@ -1,4 +1,4 @@
-version: "0.19"
+version: "0.20"
 engine: julia
 
 project:
@@ -155,4 +155,5 @@ format:
 
 execute:
   error: false
+#  freeze: false
   freeze: auto
diff --git a/quarto/alternatives/SciML.qmd b/quarto/alternatives/SciML.qmd
index 9bb6a670..38445075 100644
--- a/quarto/alternatives/SciML.qmd
+++ b/quarto/alternatives/SciML.qmd
@@ -173,7 +173,7 @@ We can solve for several parameters at once, by using a matching number of initi
 
 ```{julia}
 ps = [1, 2, 3, 4]
-u0 = @SVector[1, 1, 1, 1]
+u0 = [1, 1, 1, 1]
 prob = NonlinearProblem(f, u0, ps)
 solve(prob, NewtonRaphson())
 ```
@@ -235,6 +235,7 @@ The extra step is to specify a "`NonlinearSystem`." It is a system, as in practi
 
 ```{julia}
 ns = NonlinearSystem([eq], [x], [α], name=:ns);
+ns = complete(ns)
 ```
 
 The `name` argument is special. The name of the object (`ns`) is assigned through `=`, but the system must also know this same name. However, the name on the left is not known when the name on the right is needed, so it is up to the user to keep them synchronized. The `@named` macro handles this behind the scenes by simply rewriting the syntax of the assignment:
@@ -242,6 +243,7 @@ The `name` argument is special. The name of the object (`ns`) is assigned throug
 
 ```{julia}
 @named ns = NonlinearSystem([eq], [x], [α]);
+ns = complete(ns)
 ```
 
 With the system defined, we can pass this to `NonlinearProblem`, as was done with a function. The parameter is specified here, and in this case is `α => 1.0`. The initial guess is `[1.0]`:
@@ -366,6 +368,7 @@ The above should be self explanatory. To put into a form to pass to `solve` we d
 
 ```{julia}
 @named sys = OptimizationSystem(Area, [x], [P]);
+sys = complete(sys)
 ```
 
 (This step is different, as before an `OptimizationFunction` was defined; we use `@named`, as above, to ensure the system has the same name as the identifier, `sys`.)
diff --git a/quarto/derivatives/alternative_implicit.qmd b/quarto/derivatives/alternative_implicit.qmd
new file mode 100644
index 00000000..5526f2fe
--- /dev/null
+++ b/quarto/derivatives/alternative_implicit.qmd
@@ -0,0 +1,78 @@
+## THIS IS FAILING
+## Appendix
+
+
+There are other packages in the `Julia` ecosystem that can plot implicit equations.
+
+
+### The ImplicitEquations package
+
+
+The `ImplicitEquations` packages can plot equations and inequalities. The use is somewhat similar to the examples above, but the object plotted is a predicate, not a function. These predicates are created with functions like `Eq` or `Lt`.
+
+
+For example, the `ImplicitPlots` manual shows this function $f(x,y) = (x^4 + y^4 - 1) \cdot (x^2 + y^2 - 2) + x^5 \cdot y$ to plot. Using `ImplicitEquations`, this equation would be plotted with:
+
+
+```{julia}
+#| hold: true
+using ImplicitEquations
+f(x,y) = (x^4 + y^4 - 1) * (x^2 + y^2 - 2) + x^5 * y
+r = Eq(f, 0)  # the equation f(x,y) = 0
+plot(r)
+```
+
+Unlike `ImplicitPlots`, inequalities may be displayed:
+
+
+```{julia}
+#| hold: true
+f(x,y) = (x^4 + y^4 - 1) * (x^2 + y^2 - 2) + x^5 * y
+r = Lt(f, 0)          # the inequality f(x,y) < 0
+plot(r; M=10, N=10)  # less blocky
+```
+
+The rendered plots look "blocky" due to the algorithm used to plot the equations. As there is no rule defining $(x,y)$ pairs to plot, a search by regions is done. A region is initially labeled undetermined. If it can be shown that for any value in the region the equation is true (equations can also be inequalities), the region is colored black. If it can be shown it will never be true, the region is dropped. If a black-and-white answer is not clear, the region is subdivided and each subregion is similarly tested. This continues until the remaining undecided regions are smaller than some threshold. Such regions comprise a boundary, and here are also colored black. Only regions are plotted - not $(x,y)$ pairs - so the results are blocky. Pass larger values of $N=M$ (with defaults of $8$) to `plot` to lower the threshold at the cost of longer computation times, as seen in the last example.
+
+
+### The IntervalConstraintProgramming package
+
+
+The `IntervalConstraintProgramming` package also can be used to graph implicit equations. For certain problem descriptions it is significantly faster and makes better graphs. The usage is slightly more involved. We show the commands, but don't run them here, as there are minor conflicts with the `CalculusWithJulia`package.
+
+
+We specify a problem using the `@constraint` macro. Using a macro allows expressions to involve free symbols, so the problem is specified in an equation-like manner:
+
+
+```{julia}
+#| eval: false
+S = @constraint x^2 + y^2 <= 2
+```
+
+The right hand side must be a number.
+
+
+The area to plot over must be specified as an `IntervalBox`, basically a pair of intervals. The interval $[a,b]$ is expressed through `a..b`:
+
+
+```{julia}
+#| eval: false
+J = -3..3
+X = IntervalArithmetic.IntervalBox(J, J)
+```
+
+The `pave` command does the heavy lifting:
+
+
+```{julia}
+#| eval: false
+region = IntervalConstraintProgramming.pave(S, X)
+```
+
+A plot can be made of either the boundary, the interior, or both.
+
+
+```{julia}
+#| eval: false
+plot(region.inner)       # plot interior; use r.boundary for boundary
+```
diff --git a/quarto/derivatives/curve_sketching.qmd b/quarto/derivatives/curve_sketching.qmd
index 6ea4207b..fd792968 100644
--- a/quarto/derivatives/curve_sketching.qmd
+++ b/quarto/derivatives/curve_sketching.qmd
@@ -9,6 +9,7 @@ This section uses the following add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 using SymPy
 using Roots
 using Polynomials # some name clash with SymPy
@@ -38,7 +39,7 @@ With these, a sketch fills in between the points/lines associated with these val
 #| echo: false
 #| cache: true
 ### {{{ sketch_sin_plot }}}
-
+gr()
 
 function sketch_sin_plot_graph(i)
     f(x) = 10*sin(pi/2*x)  # [0,4]
@@ -84,7 +85,7 @@ end
 
 imgfile = tempname() * ".gif"
 gif(anim, imgfile, fps = 1)
-
+plotly()
 ImageFile(imgfile, caption)
 ```
 
diff --git a/quarto/derivatives/derivatives.qmd b/quarto/derivatives/derivatives.qmd
index d1194e10..1160ffd0 100644
--- a/quarto/derivatives/derivatives.qmd
+++ b/quarto/derivatives/derivatives.qmd
@@ -9,6 +9,7 @@ This section uses these add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 using SymPy
 ```
 
@@ -218,6 +219,7 @@ The slope of the secant line represents the average rate of change over a given
 #| hold: true
 #| echo: false
 #| cache: true
+gr()
 function secant_line_tangent_line_graph(n)
     f(x) = sin(x)
     c = pi/3
@@ -251,6 +253,7 @@ end
 imgfile = tempname() * ".gif"
 gif(anim, imgfile, fps = 1)
 
+plotly()
 ImageFile(imgfile, caption)
 ```
 
@@ -268,6 +271,7 @@ We will define the tangent line at $(c, f(c))$ to be the line through the point
 #| hold: true
 #| echo: false
 #| cache: true
+gr()
 function line_approx_fn_graph(n)
     f(x) = sin(x)
     c = pi/3
@@ -296,6 +300,7 @@ end
 imgfile = tempname() * ".gif"
 gif(anim, imgfile, fps = 1)
 
+plotly()
 ImageFile(imgfile, caption)
 ```
 
diff --git a/quarto/derivatives/first_second_derivatives.qmd b/quarto/derivatives/first_second_derivatives.qmd
index f81b5d28..9f8c749f 100644
--- a/quarto/derivatives/first_second_derivatives.qmd
+++ b/quarto/derivatives/first_second_derivatives.qmd
@@ -9,6 +9,7 @@ This section uses these add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 using SymPy
 using Roots
 ```
diff --git a/quarto/derivatives/implicit_differentiation.qmd b/quarto/derivatives/implicit_differentiation.qmd
index c565b4e0..94061878 100644
--- a/quarto/derivatives/implicit_differentiation.qmd
+++ b/quarto/derivatives/implicit_differentiation.qmd
@@ -9,6 +9,7 @@ This section uses these add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 using Roots
 using SymPy
 ```
@@ -993,81 +994,3 @@ choices = ["concave up", "concave down", "both concave up and down"]
 answ = 3
 radioq(choices, answ, keep_order=true)
 ```
-
-## Appendix
-
-
-There are other packages in the `Julia` ecosystem that can plot implicit equations.
-
-
-### The ImplicitEquations package
-
-
-The `ImplicitEquations` packages can plot equations and inequalities. The use is somewhat similar to the examples above, but the object plotted is a predicate, not a function. These predicates are created with functions like `Eq` or `Lt`.
-
-
-For example, the `ImplicitPlots` manual shows this function $f(x,y) = (x^4 + y^4 - 1) \cdot (x^2 + y^2 - 2) + x^5 \cdot y$ to plot. Using `ImplicitEquations`, this equation would be plotted with:
-
-
-```{julia}
-#| hold: true
-using ImplicitEquations
-f(x,y) = (x^4 + y^4 - 1) * (x^2 + y^2 - 2) + x^5 * y
-r = Eq(f, 0)  # the equation f(x,y) = 0
-plot(r)
-```
-
-Unlike `ImplicitPlots`, inequalities may be displayed:
-
-
-```{julia}
-#| hold: true
-f(x,y) = (x^4 + y^4 - 1) * (x^2 + y^2 - 2) + x^5 * y
-r = Lt(f, 0)          # the inequality f(x,y) < 0
-plot(r; M=10, N=10)  # less blocky
-```
-
-The rendered plots look "blocky" due to the algorithm used to plot the equations. As there is no rule defining $(x,y)$ pairs to plot, a search by regions is done. A region is initially labeled undetermined. If it can be shown that for any value in the region the equation is true (equations can also be inequalities), the region is colored black. If it can be shown it will never be true, the region is dropped. If a black-and-white answer is not clear, the region is subdivided and each subregion is similarly tested. This continues until the remaining undecided regions are smaller than some threshold. Such regions comprise a boundary, and here are also colored black. Only regions are plotted - not $(x,y)$ pairs - so the results are blocky. Pass larger values of $N=M$ (with defaults of $8$) to `plot` to lower the threshold at the cost of longer computation times, as seen in the last example.
-
-
-### The IntervalConstraintProgramming package
-
-
-The `IntervalConstraintProgramming` package also can be used to graph implicit equations. For certain problem descriptions it is significantly faster and makes better graphs. The usage is slightly more involved. We show the commands, but don't run them here, as there are minor conflicts with the `CalculusWithJulia`package.
-
-
-We specify a problem using the `@constraint` macro. Using a macro allows expressions to involve free symbols, so the problem is specified in an equation-like manner:
-
-
-```{julia}
-#| eval: false
-S = @constraint x^2 + y^2 <= 2
-```
-
-The right hand side must be a number.
-
-
-The area to plot over must be specified as an `IntervalBox`, basically a pair of intervals. The interval $[a,b]$ is expressed through `a..b`:
-
-
-```{julia}
-#| eval: false
-J = -3..3
-X = IntervalArithmetic.IntervalBox(J, J)
-```
-
-The `pave` command does the heavy lifting:
-
-
-```{julia}
-#| eval: false
-region = IntervalConstraintProgramming.pave(S, X)
-```
-
-A plot can be made of either the boundary, the interior, or both.
-
-
-```{julia}
-#| eval: false
-plot(region.inner)       # plot interior; use r.boundary for boundary
-```
diff --git a/quarto/derivatives/lhospitals_rule.qmd b/quarto/derivatives/lhospitals_rule.qmd
index 5baacfff..829a62b8 100644
--- a/quarto/derivatives/lhospitals_rule.qmd
+++ b/quarto/derivatives/lhospitals_rule.qmd
@@ -9,8 +9,8 @@ This section uses these add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 using SymPy
-
 ```
 
 
@@ -245,6 +245,7 @@ A first proof of L'Hospital's rule takes advantage of Cauchy's [generalization](
 #| echo: false
 #| cache: true
 using Roots
+gr()
 let
 ## {{{lhopitals_picture}}}
 
@@ -293,7 +294,7 @@ imgfile = tempname() * ".gif"
 gif(anim, imgfile, fps = 1)
 
 
-plotly()
+    plotly()
     ImageFile(imgfile, caption)
 end
 ```
diff --git a/quarto/derivatives/linearization.qmd b/quarto/derivatives/linearization.qmd
index 2290342a..1ab5e630 100644
--- a/quarto/derivatives/linearization.qmd
+++ b/quarto/derivatives/linearization.qmd
@@ -9,6 +9,7 @@ This section uses these add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 using SymPy
 using TaylorSeries
 using DualNumbers
diff --git a/quarto/derivatives/mean_value_theorem.qmd b/quarto/derivatives/mean_value_theorem.qmd
index 21ae30a2..318f855c 100644
--- a/quarto/derivatives/mean_value_theorem.qmd
+++ b/quarto/derivatives/mean_value_theorem.qmd
@@ -9,8 +9,8 @@ This section uses these add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 using Roots
-
 ```
 
 ```{julia}
@@ -288,13 +288,19 @@ j(x) = exp(x) * x * (x-1)
 find_zeros(j', 0, 1)
 ```
 
-This graph illustrates the lone value for $c$ for this problem
+The following graph illustrates the lone value for $c$ in $[a,b]$ for
+this problem:
 
 
 ```{julia}
 #| echo: false
 x0 = find_zero(j', (0, 1))
-plot([j, x->j(x0) + 0*(x-x0)], 0, 1)
+j₀ = j(x0)
+plot([j, x->j₀ + 0*(x-x0)], 0, 1; legend=false)
+scatter!([0,x0,1], [j₀, j₀, j₀])
+annotate!([(0,j₀,text("a", :bottom)),
+          (x0, j₀, text("c", :bottom)),
+          (1, j₀, text("b", :bottom))])
 ```
 
 ## The mean value theorem
@@ -327,10 +333,15 @@ cps = find_zeros(x -> f'(x) - m, a, b)
 p = plot(f, a-1, b+1,         linewidth=3, legend=false)
 plot!(x -> f(a) + m*(x-a),  a-1, b+1,   linewidth=3, color=:orange)
 scatter!([a,b], [f(a), f(b)])
+annotate!([(a, f(a), text("a", :bottom)),
+           (b, f(b), text("b", :bottom))])
 
 for cp in cps
-  plot!(x -> f(cp) + f'(cp)*(x-cp), a-1, b+1, color=:red)
+    plot!(x -> f(cp) + f'(cp)*(x-cp), a-1, b+1, color=:red)
 end
+scatter!(cps, f.(cps))
+subsscripts = collect("₀₁₂₃₄₅₆₇₈₉")
+annotate!([(cp, f(cp), text("c"*subsscripts[i], :bottom)) for (i,cp) ∈ enumerate(cps)])
 p
 ```
 
@@ -454,7 +465,7 @@ The Cauchy mean value theorem can be visualized in terms of a tangent line and a
 #| echo: false
 #| cache: true
 ### {{{parametric_fns}}}
-
+gr()
 
 
 function parametric_fns_graph(n)
@@ -488,7 +499,7 @@ end
 
 imgfile = tempname() * ".gif"
 gif(anim, imgfile, fps = 1)
-
+plotly()
 ImageFile(imgfile, caption)
 ```
 
diff --git a/quarto/derivatives/more_zeros.qmd b/quarto/derivatives/more_zeros.qmd
index 3932ed4f..9995045b 100644
--- a/quarto/derivatives/more_zeros.qmd
+++ b/quarto/derivatives/more_zeros.qmd
@@ -9,7 +9,7 @@ This section uses these add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
-using ImplicitEquations
+plotly()
 using Roots
 using SymPy
 ```
@@ -23,7 +23,7 @@ Newton's method is not the only algorithm of its kind for identifying zeros of a
 ## The `find_zero(f, x0)` function
 
 
-The  function `find_zero` from the `Roots` packages provides several different algorithms for finding a zero of a function, including some a derivative-free algorithms for finding zeros  when started with an initial guess. The default method is similar to Newton's method in that only a good initial guess is needed. However, the algorithm, while possibly slower in terms of function evaluations and steps,  is engineered to be a bit more robust to the choice of initial estimate than Newton's method. (If it finds a bracket, it will use a bisection algorithm which is guaranteed to converge, but can be slower to do so.) Here we see how to call the function:
+The  function `find_zero` from the `Roots` packages provides several different algorithms for finding a zero of a function, including some derivative-free algorithms for finding zeros when started with a nearby initial guess. The default method is similar to Newton's method in that only a good, initial guess is needed. However, the algorithm, while possibly slower in terms of function evaluations and steps,  is engineered to be a bit more robust to the choice of initial estimate than Newton's method. (If it finds a bracket, it will use a bisection algorithm which is guaranteed to converge, but can be slower to do so.) Here we see how to call the function:
 
 
 ```{julia}
@@ -279,10 +279,14 @@ We can see it in action on the sine function.   Here we pass in $\lambda$, but i
 chandrapatla(sin, 3, 4,  λ3, verbose=true)
 ```
 
+```{julia}
+#| output: false
+
+#=
 The condition `Φ^2 < ξ < 1 - (1-Φ)^2` can be visualized. Assume `a,b=0,1`, `fa,fb=-1/2,1`, Then `c < a < b`, and `fc` has the same sign as `fa`, but what values of `fc` will satisfy the inequality?
 
 
-```{julia}
+XX```{julia}
 ξ(c,fc) = (a-b)/(c-b)
 Φ(c,fc) = (fa-fb)/(fc-fb)
 Φl(c,fc) = Φ(c,fc)^2
@@ -291,7 +295,7 @@ a,b = 0, 1
 fa,fb = -1/2, 1
 region = Lt(Φl, ξ) & Lt(ξ,Φr)
 plot(region, xlims=(-2,a), ylims=(-3,0))
-```
+XX```
 
 When `(c,fc)` is in the shaded area, the inverse quadratic step is chosen. We can see that `fc < fa` is needed.
 
@@ -299,14 +303,16 @@ When `(c,fc)` is in the shaded area, the inverse quadratic step is chosen. We ca
 For these values, this area is within the area where a implicit quadratic step will result in a value between `a` and `b`:
 
 
-```{julia}
+XX```{julia}
 l(c,fc) = λ3(fa,fb,fc,a,b,c)
 region₃ = ImplicitEquations.Lt(l,b) & ImplicitEquations.Gt(l,a)
 plot(region₃, xlims=(-2,0), ylims=(-3,0))
-```
+XX```
 
 There are values in the parameter space where this does not occur.
-
+=#
+nothing
+```
 
 ## Tolerances
 
diff --git a/quarto/derivatives/newtons_method.qmd b/quarto/derivatives/newtons_method.qmd
index bba8b108..98631cd8 100644
--- a/quarto/derivatives/newtons_method.qmd
+++ b/quarto/derivatives/newtons_method.qmd
@@ -9,6 +9,7 @@ This section uses these add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 using SymPy
 using Roots
 ```
@@ -20,7 +21,7 @@ using Roots
 The Babylonian method is an algorithm to find an approximate value for $\sqrt{k}$. It was described by the first-century Greek mathematician Hero of [Alexandria](http://en.wikipedia.org/wiki/Babylonian_method).
 
 
-The method starts with some initial guess, called $x_0$. It then applies a formula to produce an improved guess. This is repeated until the improved guess is accurate enough or it is clear the algorithm fails to work.
+The method starts with some initial guess, called $x_0$. This is usually some nearby value to the answer. The method then applies a formula to produce an improved guess. This is repeated until the improved guess is accurate enough or it is clear the algorithm fails to work.
 
 
 For the Babylonian method, the next guess, $x_{i+1}$, is derived from the current guess, $x_i$. In mathematical notation, this is the updating step:
@@ -250,7 +251,7 @@ nothing
 #| echo: false
 #| cache: true
 ### {{{newtons_method_example}}}
-
+gr()
 caption = """
 
 Illustration of Newton's Method converging to a zero of a function.
@@ -266,7 +267,7 @@ end
 
 imgfile = tempname() * ".gif"
 gif(anim, imgfile, fps = 1)
-
+plotly()
 ImageFile(imgfile, caption)
 ```
 
@@ -730,6 +731,7 @@ What can go wrong when one of these isn't the case is illustrated next:
 #| echo: false
 #| cache: true
 ### {{{newtons_method_poor_x0}}}
+gr()
 caption = """
 
 Illustration of Newton's Method converging to a zero of a function,
@@ -748,7 +750,7 @@ end
 
 imgfile = tempname() * ".gif"
 gif(anim, imgfile, fps = 2)
-
+plotly()
 ImageFile(imgfile, caption)
 ```
 
@@ -757,6 +759,7 @@ ImageFile(imgfile, caption)
 #| echo: false
 #| cache: true
 # {{{newtons_method_flat}}}
+gr()
 caption = L"""
 
 Illustration of Newton's method failing to converge as for some $x_i$,
@@ -777,7 +780,7 @@ anim = @animate for i=1:n
 end
 imgfile = tempname() * ".gif"
 gif(anim, imgfile, fps = 1)
-
+plotly()
 ImageFile(imgfile, caption)
 ```
 
@@ -789,7 +792,7 @@ ImageFile(imgfile, caption)
 #| echo: false
 #| cache: true
 # {{{newtons_method_cycle}}}
-
+gr()
 fn, a, b, c, = x -> abs(x)^(0.49),  -2, 2, 1.0
 caption = L"""
 
@@ -808,6 +811,7 @@ end
 imgfile = tempname() * ".gif"
 gif(anim, imgfile, fps = 2)
 
+plotly()
 ImageFile(imgfile, caption)
 ```
 
@@ -819,7 +823,7 @@ ImageFile(imgfile, caption)
 #| echo: false
 #| cache: true
 # {{{newtons_method_wilkinson}}}
-
+gr()
 caption = L"""
 
 The function $f(x) = x^{20} - 1$ has two bad behaviours for Newton's
@@ -845,7 +849,7 @@ anim = @animate for i=1:n
 end
 imgfile = tempname() * ".gif"
 gif(anim, imgfile, fps = 1)
-
+plotly()
 ImageFile(imgfile, caption)
 ```
 
diff --git a/quarto/derivatives/numeric_derivatives.qmd b/quarto/derivatives/numeric_derivatives.qmd
index 6e8b76df..86069782 100644
--- a/quarto/derivatives/numeric_derivatives.qmd
+++ b/quarto/derivatives/numeric_derivatives.qmd
@@ -9,6 +9,7 @@ This section uses these add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 using ForwardDiff
 using SymPy
 using Roots
diff --git a/quarto/derivatives/optimization.qmd b/quarto/derivatives/optimization.qmd
index fc4ce2af..11aed8fb 100644
--- a/quarto/derivatives/optimization.qmd
+++ b/quarto/derivatives/optimization.qmd
@@ -9,6 +9,7 @@ This section uses these add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 using Roots
 using SymPy
 ```
@@ -60,7 +61,7 @@ We began with the question of which rectangles of perimeter $20$ have the larges
 #| cache: true
 ### {{{perimeter_area_graphic}}}
 
-
+gr()
 function perimeter_area_graphic_graph(n)
     h = 1 + 2n
     w = 10-h
@@ -83,7 +84,7 @@ end
 
 imgfile = tempname() * ".gif"
 gif(anim, imgfile, fps = 1)
-
+plotly()
 ImageFile(imgfile, caption)
 ```
 
diff --git a/quarto/derivatives/related_rates.qmd b/quarto/derivatives/related_rates.qmd
index c6f950f1..f6fd1690 100644
--- a/quarto/derivatives/related_rates.qmd
+++ b/quarto/derivatives/related_rates.qmd
@@ -9,6 +9,7 @@ This section uses these add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 using Roots
 using SymPy
 ```
@@ -36,6 +37,7 @@ The following is a typical "book" problem:
 #| cache: true
 ### {{{growing_rects}}}
 ## Secant line approaches tangent line...
+gr()
 function growing_rects_graph(n)
     w = (t) -> 2 + 4t
     h = (t) -> 3/2 * w(t)
@@ -67,7 +69,7 @@ end
 
 imgfile = tempname() * ".gif"
 gif(anim, imgfile, fps = 1)
-
+plotly()
 ImageFile(imgfile, caption)
 ```
 
@@ -400,6 +402,7 @@ $$
 #| cache: true
 ###{{{baseball_been_berry_good}}}
 ## Secant line approaches tangent line...
+gr()
 function baseball_been_berry_good_graph(n)
 
     v0 = 15
@@ -438,7 +441,7 @@ end
 
 imgfile = tempname() * ".gif"
 gif(anim, imgfile, fps = 1)
-
+plotly()
 ImageFile(imgfile, caption)
 ```
 
diff --git a/quarto/derivatives/taylor_series_polynomials.qmd b/quarto/derivatives/taylor_series_polynomials.qmd
index 539d2d51..45fefa3e 100644
--- a/quarto/derivatives/taylor_series_polynomials.qmd
+++ b/quarto/derivatives/taylor_series_polynomials.qmd
@@ -9,6 +9,7 @@ This section uses these add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 using SymPy
 using Unitful
 ```
@@ -37,6 +38,7 @@ We will see in this section how the Taylor polynomial answers these questions, a
 #| echo: false
 #| cache: true
 ###{{{taylor_animation}}}
+gr()
 taylor(f, x, c, n) = series(f, x, c, n+1).removeO()
 function make_taylor_plot(u, a, b, k)
     k = 2k
@@ -67,7 +69,7 @@ caption = L"""
 Illustration of the Taylor polynomial of degree $k$, $T_k(x)$, at $c=0$ and its graph overlaid on that of the function $1 - \cos(x)$.
 
 """
-
+plotly()
 ImageFile(imgfile, caption)
 ```
 
diff --git a/quarto/differentiable_vector_calculus/plots_plotting.qmd b/quarto/differentiable_vector_calculus/plots_plotting.qmd
index 0d0c39ed..fe90ffb5 100644
--- a/quarto/differentiable_vector_calculus/plots_plotting.qmd
+++ b/quarto/differentiable_vector_calculus/plots_plotting.qmd
@@ -9,6 +9,7 @@ This section uses these add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 import Contour: contours, levels, level, lines, coordinates
 using LinearAlgebra
 using ForwardDiff
@@ -359,7 +360,6 @@ We can parameterize the sphere by plotting values for $x$, $y$, and $z$ produced
 
 
 ```{julia}
-#| eval: false
 thetas = range(0, stop=pi,   length=50)
 phis   = range(0, stop=pi/2, length=50)
 
@@ -379,8 +379,6 @@ For convenience, the `plot_parametric` function from `CalculusWithJulia` can pro
 
 
 ```{julia}
-#| eval: false
-#| hold: true
 F(theta, phi) = [X(1, theta, phi), Y(1, theta, phi), Z(1, theta, phi)]
 plot_parametric(0..pi, 0..pi/2, F)
 ```
@@ -396,13 +394,12 @@ $$
 where $t \in [0,1]$ This gives a vector-valued function
 
 $$
-F(u, t) = t * (r(u) - V) + V, \quad \alpha \leq u \leq \beta, 0 \leq t \leq 1.
+F(u, t) = t \cdot (r(u) - V) + V, \quad \alpha \leq u \leq \beta, 0 \leq t \leq 1.
 $$
 
 To illustrate, we have:
 
 ```{julia}
-#| eval: false
 # cf. https://discourse.julialang.org/t/general-plotting-code-for-cone-in-3d-with-glmakie-or-plots/92104/3
 
 basecurve(u) = [cos(u), sin(u) + sin(u/2), 0]
@@ -421,7 +418,6 @@ To use it, we see what happens when a sphere if rendered:
 
 
 ```{julia}
-#| eval: false
 #| hold: true
 f(x,y,z) = x^2 + y^2 + z^2 - 25
 CalculusWithJulia.plot_implicit_surface(f)
@@ -431,7 +427,6 @@ This figure comes from a February 14, 2019 article in the [New York Times](https
 
 
 ```{julia}
-#| eval: false
 #| hold: true
 a,b = 1,3
 f(x,y,z) = (x^2+((1+b)*y)^2+z^2-1)^3-x^2*z^3-a*y^2*z^3
diff --git a/quarto/differentiable_vector_calculus/polar_coordinates.qmd b/quarto/differentiable_vector_calculus/polar_coordinates.qmd
index 48c384cc..d1bc568d 100644
--- a/quarto/differentiable_vector_calculus/polar_coordinates.qmd
+++ b/quarto/differentiable_vector_calculus/polar_coordinates.qmd
@@ -9,6 +9,7 @@ This section uses these add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 using SymPy
 using Roots
 using QuadGK
@@ -29,6 +30,7 @@ Polar coordinates parameterize the plane though an angle $\theta$ made from the
 ```{julia}
 #| hold: true
 #| echo: false
+gr()
 theta = pi/6
 rr = 1
 
@@ -48,7 +50,9 @@ annotate!([
         (cos(theta), sin(theta)+.05, L"(x,y)"),
         (cos(theta),-.05, L"x"),
         (-.05, sin(theta),L"y")
-        ])
+])
+plotly()
+p
 ```
 
 To recover the Cartesian coordinates from the pair $(r,\theta)$, we have these formulas from [right](http://en.wikipedia.org/wiki/Polar_coordinate_system#Converting_between_polar_and_Cartesian_coordinates) triangle geometry:
diff --git a/quarto/differentiable_vector_calculus/scalar_functions.qmd b/quarto/differentiable_vector_calculus/scalar_functions.qmd
index bccd006d..91e55df2 100644
--- a/quarto/differentiable_vector_calculus/scalar_functions.qmd
+++ b/quarto/differentiable_vector_calculus/scalar_functions.qmd
@@ -9,6 +9,7 @@ This section uses these add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+gr()
 using ForwardDiff
 using SymPy
 using Roots
@@ -208,7 +209,7 @@ An alternate to `surface` is `wireframe` – which may not use shading in all ba
 ```{julia}
 #| hold: true
 xs = ys = range(-2,2, length=10) # downsample to see the frame
-wireframe(xs, ys, 𝒇)   # gr() or pyplot() wireplots render better than plotly()
+wireframe(xs, ys, 𝒇)
 ```
 
 ##### Example
diff --git a/quarto/differentiable_vector_calculus/scalar_functions_applications.qmd b/quarto/differentiable_vector_calculus/scalar_functions_applications.qmd
index ed49f111..2c6ac567 100644
--- a/quarto/differentiable_vector_calculus/scalar_functions_applications.qmd
+++ b/quarto/differentiable_vector_calculus/scalar_functions_applications.qmd
@@ -9,6 +9,7 @@ This section uses these add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 using SymPy
 using Roots
 ```
diff --git a/quarto/differentiable_vector_calculus/vector_fields.qmd b/quarto/differentiable_vector_calculus/vector_fields.qmd
index c344ec79..c4daf38a 100644
--- a/quarto/differentiable_vector_calculus/vector_fields.qmd
+++ b/quarto/differentiable_vector_calculus/vector_fields.qmd
@@ -9,6 +9,7 @@ This section uses these add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 using SymPy
 using ForwardDiff
 using LinearAlgebra
diff --git a/quarto/differentiable_vector_calculus/vector_valued_functions.qmd b/quarto/differentiable_vector_calculus/vector_valued_functions.qmd
index 2ea85c41..a382a328 100644
--- a/quarto/differentiable_vector_calculus/vector_valued_functions.qmd
+++ b/quarto/differentiable_vector_calculus/vector_valued_functions.qmd
@@ -9,6 +9,7 @@ This section uses these add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 using SymPy
 using Roots
 using LinearAlgebra
@@ -322,11 +323,13 @@ In 1935 [Marcel Duchamp](https://arthur.io/art/marcel-duchamp/rotorelief-no-10-c
 
 
 ```{julia}
+#| echo: false
 let
 # https://exploratorium.tumblr.com/post/33140874462/marcel-duchamp-rotoreliefs-duchamp-recognized
 
 # coordinates and colors selected by gimp from
-# https://arthur.io/art/marcel-duchamp/rotorelief-no-10-cage-modele-depose-verso
+    # https://arthur.io/art/marcel-duchamp/rotorelief-no-10-cage-modele-depose-verso
+    gr()
     circs = [466 548 513 505 556 554 # x₁,y₁,x₂,y₂,x₂,y₃
              414 549 511 455 595 549
              365 545 507 408 635 548
@@ -354,7 +357,7 @@ let
     function solve_i(i)
         eqs = [(p[1] - h)^2 + (p[2]-k)^2 ~ r^2 for
 		       p ∈ (circs[i,1:2], circs[i,3:4], circs[i,5:6])]
-        d = solve(eqs)[1]
+        d = solve(tuple(eqs...))[1]
         (x=float(d[h]), y=float(d[k]), r=float(d[r]))
     end
     c₀, cs... = solve_i.(16:-1:1) # c₀ is centered
@@ -391,6 +394,9 @@ let
 
     fname = tempname() * ".gif"
     gif(anim, fname, fps = 40)
+
+    plotly()
+    ImageFile(fname, "Duchamp rotorelief")
 end
 ```
 
@@ -1268,6 +1274,10 @@ For the torsion to be defined, the cross product $\vec{r}' \times \vec{r}''$ mus
 
 ##### Example: Tubular surface
 
+```{julia}
+#| echo: false
+gr();
+```
 
 This last example comes from a collection of several [examples](https://github.com/empet/3D-Viz-with-PlotlyJS.jl/blob/main/5-Tubular-surface.ipynb) provided by Discourse user `@empet` to illustrate `plotlyjs`. We adopt it to `Plots` with some minor changes below.
 
@@ -1376,6 +1386,11 @@ ts = range(t_span..., length=1001)
 surface(unzip(r_0.(ts, θs'))...)
 ```
 
+```{julia}
+#| echo: false
+plotly();
+```
+
 ## Arc length
 
 
@@ -1586,6 +1601,7 @@ Here, when $b$ gets large, the curve looks more and more "straight" and the tors
 ```{julia}
 #| echo: false
 let
+    gr()
 a = 1
 F(t) = [cos(pi/2 - t), 2sin(pi/2-t)]
 p = (a, F)
@@ -1613,7 +1629,7 @@ using LaTeXStrings
 annotate!([(-.5,1.5,L"k"),
 (.775,1.55,L"\kappa"),
 (.85, 1.3, L"\alpha")])
-
+plotly()
 plt
 end
 ```
diff --git a/quarto/integral_vector_calculus/div_grad_curl.qmd b/quarto/integral_vector_calculus/div_grad_curl.qmd
index d3ff059c..f043b8c1 100644
--- a/quarto/integral_vector_calculus/div_grad_curl.qmd
+++ b/quarto/integral_vector_calculus/div_grad_curl.qmd
@@ -9,6 +9,7 @@ This section uses these add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 using SymPy
 ```
 
diff --git a/quarto/integral_vector_calculus/double_triple_integrals.qmd b/quarto/integral_vector_calculus/double_triple_integrals.qmd
index fcab403f..8f828899 100644
--- a/quarto/integral_vector_calculus/double_triple_integrals.qmd
+++ b/quarto/integral_vector_calculus/double_triple_integrals.qmd
@@ -9,6 +9,7 @@ This section uses these add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 using QuadGK
 using SymPy
 using HCubature
diff --git a/quarto/integral_vector_calculus/line_integrals.qmd b/quarto/integral_vector_calculus/line_integrals.qmd
index 56e02540..6e246749 100644
--- a/quarto/integral_vector_calculus/line_integrals.qmd
+++ b/quarto/integral_vector_calculus/line_integrals.qmd
@@ -9,6 +9,7 @@ This section uses these add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 using QuadGK
 using SymPy
 using HCubature
diff --git a/quarto/integral_vector_calculus/review.qmd b/quarto/integral_vector_calculus/review.qmd
index 36e0b8a3..812a06e6 100644
--- a/quarto/integral_vector_calculus/review.qmd
+++ b/quarto/integral_vector_calculus/review.qmd
@@ -8,6 +8,7 @@
 #| results: "hidden"
 using CalculusWithJulia
 using Plots
+plotly()
 nothing
 ```
 
diff --git a/quarto/integral_vector_calculus/stokes_theorem.qmd b/quarto/integral_vector_calculus/stokes_theorem.qmd
index 3f7de1c6..afd2bbaa 100644
--- a/quarto/integral_vector_calculus/stokes_theorem.qmd
+++ b/quarto/integral_vector_calculus/stokes_theorem.qmd
@@ -9,6 +9,7 @@ This section uses these add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 using QuadGK
 using SymPy
 ```
diff --git a/quarto/integrals/arc_length.qmd b/quarto/integrals/arc_length.qmd
index d0315612..10d164cb 100644
--- a/quarto/integrals/arc_length.qmd
+++ b/quarto/integrals/arc_length.qmd
@@ -9,6 +9,7 @@ This section uses these add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 using SymPy
 using QuadGK
 using Roots
@@ -69,6 +70,7 @@ To see why, any partition of the interval $[a,b]$ by $a = t_0 < t_1 < \cdots < t
 #| hold: false
 #| echo: false
 ## {{{arclength_graph}}}
+gr()
 function make_arclength_graph(n)
 
     ns = [10,15,20, 30, 50]
@@ -100,6 +102,7 @@ The arc length of the parametric curve can be approximated using straight line s
 
 """
 
+plotly()
 ImageFile(imgfile, caption)
 ```
 
diff --git a/quarto/integrals/area.qmd b/quarto/integrals/area.qmd
index 1e5326f6..c6546522 100644
--- a/quarto/integrals/area.qmd
+++ b/quarto/integrals/area.qmd
@@ -9,6 +9,7 @@ This section uses these add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 using QuadGK
 using Roots
 ```
@@ -75,7 +76,7 @@ In a previous section, we saw this animation:
 #| hold: true
 #| echo: false
 ## {{{archimedes_parabola}}}
-
+gr()
 
 f(x) = x^2
 colors = [:black, :blue, :orange, :red, :green, :orange, :purple]
@@ -126,7 +127,7 @@ caption = L"""
 The first triangle has area $1/2$, the second has area $1/8$, then $2$ have area $(1/8)^2$, $4$ have area $(1/8)^3$, ...
 With some algebra, the total area then should be $1/2 \cdot (1 + (1/4) + (1/4)^2 + \cdots) = 2/3$.
 """
-
+plotly()
 ImageFile(imgfile, caption)
 ```
 
@@ -188,6 +189,7 @@ Clearly for a given partition and choice of $c_i$, the above can be computed. Ea
 ```{julia}
 #| hold: true
 #| echo: false
+gr()
 rectangle(x, y, w, h) = Shape(x .+ [0,w,w,0], y .+ [0,0,h,h])
 function ₙ(j)
 	a = ("₋","","","₀","₁","₂","₃","₄","₅","₆","₇","₈","₉")
@@ -218,6 +220,7 @@ gif(anim, imgfile, fps = 1)
 
 caption = "Illustration of left Riemann sum for increasing ``n`` values"
 
+plotly()
 ImageFile(imgfile, caption)
 ```
 
diff --git a/quarto/integrals/area_between_curves.qmd b/quarto/integrals/area_between_curves.qmd
index 559b6f2a..bedb2b4b 100644
--- a/quarto/integrals/area_between_curves.qmd
+++ b/quarto/integrals/area_between_curves.qmd
@@ -9,6 +9,7 @@ This section uses these add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 using Roots
 using QuadGK
 using SymPy
diff --git a/quarto/integrals/center_of_mass.qmd b/quarto/integrals/center_of_mass.qmd
index 3b30c316..5c0f1cbc 100644
--- a/quarto/integrals/center_of_mass.qmd
+++ b/quarto/integrals/center_of_mass.qmd
@@ -9,6 +9,7 @@ This section uses these add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 using Roots
 using QuadGK
 using SymPy
diff --git a/quarto/integrals/ftc.qmd b/quarto/integrals/ftc.qmd
index bc9c5a64..f1e38b48 100644
--- a/quarto/integrals/ftc.qmd
+++ b/quarto/integrals/ftc.qmd
@@ -9,6 +9,7 @@ This section uses these add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 using SymPy
 using Roots
 using QuadGK
@@ -434,7 +435,7 @@ The value of $a$ does not matter, as long as the integral is defined.
 #| hold: true
 #| echo: false
 ##{{{ftc_graph}}}
-
+gr()
 function make_ftc_graph(n)
     a, b = 2, 3
     ts = range(0, stop=b, length=50)
@@ -471,7 +472,7 @@ end
 
 imgfile = tempname() * ".gif"
 gif(anim, imgfile, fps = 1)
-
+plotly()
 ImageFile(imgfile, caption)
 ```
 
diff --git a/quarto/integrals/improper_integrals.qmd b/quarto/integrals/improper_integrals.qmd
index 0f8f6891..b0927021 100644
--- a/quarto/integrals/improper_integrals.qmd
+++ b/quarto/integrals/improper_integrals.qmd
@@ -9,6 +9,7 @@ This section uses these add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 using SymPy
 using QuadGK
 ```
@@ -27,7 +28,7 @@ To define integrals  with either functions having singularities or infinite doma
 #| hold: true
 #| echo: false
 ### {{{sqrt_graph}}}
-
+gr()
 function make_sqrt_x_graph(n)
 
     b = 1
@@ -58,7 +59,7 @@ end
 
 imgfile = tempname() * ".gif"
 gif(anim, imgfile, fps = 1)
-
+plotly()
 ImageFile(imgfile, caption)
 ```
 
diff --git a/quarto/integrals/integration_by_parts.qmd b/quarto/integrals/integration_by_parts.qmd
index 1c5ebdc5..7e99c021 100644
--- a/quarto/integrals/integration_by_parts.qmd
+++ b/quarto/integrals/integration_by_parts.qmd
@@ -9,6 +9,7 @@ This section uses these add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 using SymPy
 ```
 
diff --git a/quarto/integrals/mean_value_theorem.qmd b/quarto/integrals/mean_value_theorem.qmd
index ef03c2ae..e28bba76 100644
--- a/quarto/integrals/mean_value_theorem.qmd
+++ b/quarto/integrals/mean_value_theorem.qmd
@@ -9,6 +9,7 @@ This section uses these add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 using QuadGK
 ```
 
diff --git a/quarto/integrals/substitution.qmd b/quarto/integrals/substitution.qmd
index 0578811e..a3dd9634 100644
--- a/quarto/integrals/substitution.qmd
+++ b/quarto/integrals/substitution.qmd
@@ -9,6 +9,7 @@ This section uses these add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 using SymPy
 
 ```
diff --git a/quarto/integrals/surface_area.qmd b/quarto/integrals/surface_area.qmd
index cec88c46..688de64f 100644
--- a/quarto/integrals/surface_area.qmd
+++ b/quarto/integrals/surface_area.qmd
@@ -10,6 +10,7 @@ This section uses these add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 using SymPy
 using QuadGK
 ```
@@ -549,7 +550,7 @@ numericq(val)
 ```{julia}
 #| hold: true
 #| echo: false
-
+gr()
 ## For **some reason** having this in the natural place messes up the plots.
 ## {{{approximate_surface_area}}}
 
@@ -582,5 +583,6 @@ Surface of revolution of $f(x) = 2 - x^2$ about the $y$ axis. The lines segments
 
 """
 
+plotly()
 ImageFile(imgfile, caption)
 ```
diff --git a/quarto/integrals/twelve-qs.qmd b/quarto/integrals/twelve-qs.qmd
index 472d423b..324432dd 100644
--- a/quarto/integrals/twelve-qs.qmd
+++ b/quarto/integrals/twelve-qs.qmd
@@ -5,6 +5,7 @@ This section uses these packages:
 ```{julia}
 using SymPy
 using Plots
+plotly()
 ```
 
 ----
diff --git a/quarto/integrals/volumes_slice.qmd b/quarto/integrals/volumes_slice.qmd
index 1156cb0d..2a891ffa 100644
--- a/quarto/integrals/volumes_slice.qmd
+++ b/quarto/integrals/volumes_slice.qmd
@@ -9,6 +9,7 @@ This section uses these add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 using QuadGK
 using Unitful, UnitfulUS
 using Roots
diff --git a/quarto/limits/continuity.qmd b/quarto/limits/continuity.qmd
index 8935cd74..6371441c 100644
--- a/quarto/limits/continuity.qmd
+++ b/quarto/limits/continuity.qmd
@@ -9,6 +9,7 @@ This section uses these add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 using SymPy
 ```
 
diff --git a/quarto/limits/intermediate_value_theorem.qmd b/quarto/limits/intermediate_value_theorem.qmd
index d48b27c3..f3b88df0 100644
--- a/quarto/limits/intermediate_value_theorem.qmd
+++ b/quarto/limits/intermediate_value_theorem.qmd
@@ -9,6 +9,7 @@ This section uses these add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 using Roots
 using SymPy
 ```
@@ -31,7 +32,7 @@ Continuity for functions is a valued property which carries implications. In thi
 #| echo: false
 #| cache: true
 ### {{{IVT}}}
-
+gr()
 function IVT_graph(n)
     f(x) = sin(pi*x) + 9x/10
     a,b = [0,3]
@@ -68,6 +69,7 @@ with $f(x)=y$.
 
 """
 
+plotly()
 ImageFile(imgfile, caption)
 ```
 
@@ -88,12 +90,20 @@ We use this fact when building a "sign chart" of a polynomial function. Between
 
 Here, we use the Bolzano theorem to give an algorithm - the *bisection method* - to locate the value $c$ under the assumption $f$ is continuous on $[a,b]$ and changes sign between $a$ and $b$.
 
+::: {.callout-note}
+#### Between
+
+The bisection method is used to find a zero, $c$, of $f(x)$ *between* two values, $a$ and $b$. The method is guaranteed to work under assumptions, the most important being the continuous function having different signs at $a$ and $b$.
+:::
+
+
 
 ```{julia}
 #| hold: true
 #| echo: false
 #| cache: true
 ## {{{bisection_graph}}}
+gr()
 function bisecting_graph(n)
     f(x) = x^2 - 2
     a,b = [0,2]
@@ -138,6 +148,7 @@ that the intermediate value theorem guaratees a zero.
 
 """
 
+plotly()
 ImageFile(imgfile, caption)
 ```
 
@@ -269,6 +280,13 @@ Find the zero numerically. The plot shows $q(5) < 0 < q(10)$, so $[5,10]$ is a b
 find_zero(q, (5, 10))
 ```
 
+::: {.callout-note}
+### Between need not be near
+
+Later, we will see more efficient algorithms to find a zero *near* a given guess. The bisection method finds a zero *between* two values of a bracketing interval. This interval need not be small. Indeed in many cases it can be infinite. For this particular problem, any interval like `(2,N)` will work as long as `N` is bigger than the zero and small enough that `q(N)` is finite *or* infinite *but* not `NaN`. (Basically, `q` must evaluate to a number with a sign. Here, the value of `q(Inf)` is `NaN` as it evaluates to the indeterminate `Inf - Inf`. But `q` is still not `NaN` for quite large numbers, such as `1e77`, as `x^4` can as big as `1e308` -- technically `floatmax(Float64)` -- and be finite.)
+
+:::
+
 ##### Example
 
 
@@ -317,6 +335,12 @@ h(x) = f(x) - g(x)
 find_zero(h, (0, 2))
 ```
 
+::: {.callout-note}
+### Solving `f(x) = g(x)` and `f(x) = c`
+
+The above shows a means to translate a given problem into one that can be solved with `find_zero`. Basically to solve either when a function is a non-zero constant or when a function is equal to some other function, the difference between the two sides is formed and turned into a function, called `h` above.
+:::
+
 ##### Example: Inverse functions
 
 If $f(x)$ is *monotonic* and *continuous* over an interval $[a,b]$ then it has an *inverse function*. That is for any $y$ between $f(a)$ and $f(b)$ we can find an $x$ satisfying $y = f(x)$ with $a \leq x \leq b$. This is due, of course, to both the intermediate value theorem (which guarantees an $x$) and monotonicity (which guarantees just one $x$).
@@ -379,11 +403,12 @@ plot!(plan2)
 We can see the intersection point is around $14$ and that if a family generates between $0$-$14$ bags of trash per month that plan $2$ would be cheaper.
 
 
-Let's get a numeric value, using a simple bracket and an anonymous function:
+Let's get a numeric value. We first, form an auxiliary function returning the difference:
 
 
 ```{julia}
-find_zero(x -> plan1(x) - plan2(x), (10, 20))
+h(x) = plan1(x) - plan2(x)
+find_zero(h, (10, 20))
 ```
 
 ##### Example, the flight of an arrow
@@ -495,17 +520,17 @@ plot!(d, 0, x1, label="windy day")
 The Bolzano theorem assumes a continuous function $f$, and when applicable, yields an algorithm to find a guaranteed zero.
 
 
-However, the algorithm itself does not know that the function is continuous or not, only that the function changes sign. As such, it can produce answers that are not "zeros" when used with discontinuous functions.
+However, the algorithm itself does not know that the function is continuous or not, only that the function changes sign. As such, it can produce useful answers that are not "zeros" when applied to discontinuous functions.
 
 
-In general a function over floating point values could be considered as a large table of mappings: each of the $2^{64}$ floating point values gets assigned a value. This is discrete mapping, there is nothing the computer sees related to continuity.
+In general a function over floating point values could be considered as a large table of mappings: each of the $2^{64}$ floating point values gets assigned a value. This is a discrete mapping, there is nothing the computer sees related to continuity.
 
 
 > The concept of continuity, if needed, must be verified by the user of the algorithm.
 
 
 
-We have seen this when plotting rational functions or functions with vertical asymptotes. The default algorithms just connect points with lines. The user must manage the discontinuity (by assigning some values `NaN`, say); the algorithms used do not.
+We have remarked on this discreteness when plotting rational functions or functions with vertical asymptotes. The default algorithms just connect points with lines. The user must manage the discontinuity (by assigning some values `NaN`, say); the algorithms used may not.
 
 
 In this particular case, the bisection algorithm can still be fruitful even when the function is not continuous, as the algorithm will yield information about crossing values of $0$, possibly at discontinuities. But the user of the algorithm must be aware that the answers are only guaranteed to be zeros of the function if the function is continuous and the algorithm did not check for that assumption.
@@ -579,7 +604,7 @@ f(x) = x^5 - x^4 + x^3 - x^2 + 1
 find_zeros(f, (-10, 10))
 ```
 
-(Here $-10$ and $10$ were arbitrarily chosen. Cauchy's method could be used to be more systematic.)
+(The choice of $-10$ and $10$ was arbitrary. Cauchy's method could be used to be more systematic.)
 
 
 ##### Example: Solving f(x) = g(x)
@@ -592,15 +617,15 @@ To proceed with `find_zeros`, we define $f(x) = e^x - x^5$, as $f(x) = 0$ precis
 
 
 ```{julia}
-f₁(x) = exp(x) - x^5
-zs = find_zeros(f₁, (-20,20))
+f(x) = exp(x) - x^5
+zs = find_zeros(f, (-20,20))
 ```
 
 The output of `find_zeros` is a vector of values. To check that each value is an approximate zero can be done with the "." (broadcast) syntax:
 
 
 ```{julia}
-f₁.(zs)
+f.(zs)
 ```
 
 (For a continuous function this should be the case that the values returned by `find_zeros` are approximate zeros. Bear in mind that if $f$ is not continuous the algorithm might find jumping points that are not zeros and may not even be in the domain of the function.)
@@ -612,15 +637,15 @@ f₁.(zs)
 The `find_zero` function in the `Roots` package is an interface to one of several methods. For now we focus on the *bracketing* methods, later we will see others. Bracketing methods, among others,  include `Roots.Bisection()`, the basic bisection method though with a different sense of "middle" than $(a+b)/2$ and used by default above; `Roots.A42()`, which will typically converge much faster than simple bisection; `Roots.Brent()` for the classic method of Brent, and `FalsePosition()` for a family of *regula falsi* methods. These can all be used by specifying the method in a call to `find_zero`.
 
 
-Alternatively, `Roots` implements the `CommonSolve` interface popularized by its use in the `DifferentialEquations.jl` ecosystem, a wildly successful area for `Julia`. The basic setup is two steps: setup a "problem", solve the problem.
+Alternatively, `Roots` implements the `CommonSolve` interface popularized by its use in the `DifferentialEquations.jl` ecosystem, a wildly successful area for `Julia`. The basic setup involves two steps: setup a "problem;" solve the problem.
 
 
-To set up a problem, we call `ZeroProblem` with the function and an initial interval, as in:
+To set up a problem we call `ZeroProblem` with the function and an initial interval, as in:
 
 
 ```{julia}
-f₅(x) = x^5 - x - 1
-prob = ZeroProblem(f₅, (1,2))
+f(x) = x^5 - x - 1
+prob = ZeroProblem(f, (1,2))
 ```
 
 Then we can "solve" this problem with `solve`. For example:
@@ -715,7 +740,7 @@ plot(x -> x * exp(-x), 0, 5)
 ##### Example
 
 
-The tangent function does not have a *guarantee* of absolute maximum or minimum over $(-\pi/2, \pi/2),$ as it is not *continuous* at the endpoints. In fact, it doesn't have either extrema - it has vertical asymptotes at each endpoint of this interval.
+The tangent function does not have a *guarantee* of an absolute maximum or an minimum over $(-\pi/2, \pi/2),$ as it is not *continuous* at the endpoints. In fact, it doesn't have either extrema - it has vertical asymptotes at each endpoint of this interval.
 
 
 ##### Example
@@ -729,7 +754,7 @@ The function $f(x) = x^{2/3}$ over the interval $[-2,2]$ has cusp at $0$. Howeve
 plot(x -> (x^2)^(1/3), -2, 2)
 ```
 
-(The use of just `x^(2/3)` would fail, can you guess why?)
+(The use of just `x^(2/3)` in the above would fail, can you guess why?)
 
 
 ##### Example
@@ -772,6 +797,60 @@ The *preimage* of an open set, $I$, is $\{a: f(a) \in I\}$. (All $a$ with an ima
 
 ## Questions
 
+###### Question
+Consider the following plot
+
+```{julia}
+plot(airy, -5, 5; xticks=-5:5)
+```
+
+There is a guaranteed zero between:
+
+```{julia}
+#| echo: false
+choices = ["``-4`` and ``-3``",
+           "``-3`` and ``-2``",
+           "``-2`` and ``-1``",
+           "``-1`` and ``-0``"]
+answer = 2
+radioq(choices, answer; keep_order=true)
+```
+
+
+
+###### Question
+Consider the following plot
+
+```{julia}
+plot(erf, -5, 5; xticks=-5:5)
+```
+
+There is a guaranteed zero to `erf(x) = 0.5` between
+
+
+```{julia}
+#| echo: false
+choices = ["``-1`` and ``0``",
+           "``0`` and ``1``",
+           "``1`` and ``2``"]
+answer = 2
+radioq(choices, answer; keep_order=true)
+```
+
+###### Question
+
+The hyperbolic tangent has a domain of all real numbers and a range between $-1$ and $1$. It is a strictly increasing function.
+
+Is the interval $(-\infty, \infty)$ a *bracketing* interval for the function `f(x) = tanh(x) - 1/2`?
+
+```{julia}
+#| echo: false
+choices = ["Technically no, as you can't evaluate the function at these points without using a limit; but practically yes. As the hyperbolic tangent function is monotonic and increases from ``-1`` to ``1`` in that interval, it must cross the line ``y=1/2``.",
+           "No, the function has a jump"]
+answer = 1
+radioq(choices, answer; keep_order=true)
+```
+
 
 ###### Question
 
diff --git a/quarto/limits/limits.qmd b/quarto/limits/limits.qmd
index dce00957..95469be6 100644
--- a/quarto/limits/limits.qmd
+++ b/quarto/limits/limits.qmd
@@ -9,6 +9,7 @@ This section uses the following add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 using Richardson  # for extrapolation
 using SymPy       # for symbolic limits
 ```
@@ -28,7 +29,7 @@ There wasn't a ready-made formula for the area of this shape, as was known for a
 #| echo: false
 #| cache: true
 ###{{{archimedes_parabola}}}
-
+gr()
 f(x) = x^2
 colors = [:black, :blue, :orange, :red, :green, :orange, :purple]
 
@@ -76,6 +77,7 @@ The first triangle has area $1/2$, the second has area $1/8$, then $2$ have area
 With some algebra, the total area then should be $1/2 \cdot (1 + (1/4) + (1/4)^2 + \cdots) = 2/3$.
 """
 
+plotly()
 ImageFile(imgfile, caption)
 ```
 
@@ -102,7 +104,7 @@ Next, we illustrate how Archimedes approximated $\pi$ – the ratio of the circu
 #| hold: true
 #| echo: false
 ## Archimedes approximation for pi
-
+gr()
 blue, green, purple, red = :royalblue, :forestgreen, :mediumorchid3, :brown3
 
 
@@ -161,7 +163,7 @@ end
 caption = L"""
 The ratio of the circumference of a circle to its diameter, $\pi$, can be approximated from above and below by computing the perimeters of the inscribed $n$-gons. Archimedes computed the perimeters for $n$ being $12$, $24$, $48$, and $96$ to determine that $3~1/7 \leq \pi \leq 3~10/71$.
 """
-
+plotly()
 ImageFile(p, caption)
 ```
 
@@ -217,6 +219,7 @@ The above bound comes from this figure, for small $x > 0$:
 ```{julia}
 #| hold: true
 #| echo: false
+gr()
 p = plot(x -> sqrt(1 - x^2), 0, 1, legend=false, aspect_ratio=:equal,
      linewidth=3, color=:black)
 θ = π/6
@@ -237,6 +240,7 @@ imgfile = tempname() * ".png"
 savefig(p, imgfile)
 caption = "Triangle ``ABD` has less area than the shaded wedge, which has less area than triangle ``ACD``. Their respective areas are ``(1/2)\\sin(\\theta)``, ``(1/2)\\theta``, and ``(1/2)\\tan(\\theta)``. The inequality used to show ``\\sin(x)/x`` is bounded below by ``\\cos(x)`` and above by ``1`` comes from a division by ``(1/2) \\sin(x)`` and taking reciprocals.
 "
+plotly()
 ImageFile(imgfile, caption)
 ```
 
@@ -1022,7 +1026,7 @@ Then the limit of $f$ must also  be $L$.
 ```{julia}
 #| hold: true
 #| echo: false
-
+gr()
 function squeeze_example(x)
 	x₀ = 0.5
 	plot(cos, 0, x₀, label="cos")
@@ -1043,6 +1047,7 @@ gif(anim, imgfile, fps = 1)
 caption = """
 As ``x`` goes to ``0``, the values of ``\\sin(x)/x`` are squeezed between ``\\cos(x)`` and ``1`` which both converge to ``1``.
 """
+plotly()
 ImageFile(imgfile, caption)
 ```
 
@@ -1089,6 +1094,7 @@ In this case, $\delta$ is easy to guess, as the function is linear and has slope
 #| echo: false
 #| cache: true
 ## {{{ limit_e_d }}}
+gr()
 function make_limit_e_d(n)
     f(x) = x^3
 
@@ -1135,7 +1141,7 @@ does not leave the top or bottom of the box (except possibly at the
 edges). In this example $\delta^3=\epsilon$.
 
 """
-
+plotly()
 ImageFile(imgfile, caption)
 ```
 
@@ -1247,6 +1253,7 @@ Select the graph for which there is no limit at $a$.
 ```{julia}
 #| hold: true
 #| echo: false
+gr()
 let
     p1 = plot(;axis=nothing, legend=false)
     title!(p1, "(a)")
@@ -1281,6 +1288,11 @@ let
 end
 ```
 
+```{julia}
+#| echo: false
+plotly();
+```
+
 ###### Question
 
 
diff --git a/quarto/limits/limits_extensions.qmd b/quarto/limits/limits_extensions.qmd
index a6a85b0c..83284006 100644
--- a/quarto/limits/limits_extensions.qmd
+++ b/quarto/limits/limits_extensions.qmd
@@ -9,6 +9,7 @@ This section uses the following add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 using SymPy
 ```
 
@@ -624,6 +625,8 @@ Select the graph for which the limit at $a$ is infinite.
 ```{julia}
 #| hold: true
 #| echo: false
+gr()
+
 p1 = plot(;axis=nothing, legend=false)
 title!(p1, "(a)")
 plot!(p1, x -> x^2, 0, 2, color=:black)
@@ -653,6 +656,7 @@ l = @layout[a b; c d]
 p = plot(p1, p2, p3, p4, layout=l)
 imgfile = tempname() * ".png"
 savefig(p, imgfile)
+plotly()
 hotspotq(imgfile, (1/2,1), (1/2,1))
 ```
 
@@ -665,6 +669,7 @@ Select the graph for which the limit at $\infty$ appears to be defined.
 ```{julia}
 #| hold: true
 #| echo: false
+gr()
 p1 = plot(;axis=nothing, legend=false)
 title!(p1, "(a)")
 plot!(p1, x -> x^2, 0, 2, color=:black)
@@ -690,6 +695,7 @@ l = @layout[a b; c d]
 p = plot(p1, p2, p3, p4, layout=l)
 imgfile = tempname() * ".png"
 savefig(p, imgfile)
+plotly()
 hotspotq(imgfile, (1/2,1), (1/2,1))
 ```
 
diff --git a/quarto/precalc/calculator.qmd b/quarto/precalc/calculator.qmd
index a0c68cea..5efb000f 100644
--- a/quarto/precalc/calculator.qmd
+++ b/quarto/precalc/calculator.qmd
@@ -762,6 +762,30 @@ numericq(val)
 
 ###### Question
 
+You are asked to cook chicken is an unfamiliar kitchen.  Your recipe says to turn the oven to 200 Celsius, but you the oven is calibrated in Fahrenheit. Which value is closest?
+
+```{julia}
+#| echo: false
+choices = [350, 375, 400, 425, 450]
+v = 9/5*200 + 32
+_, answer = findmin(abs.(choices .- v))
+radioq(choices, answer; keep_order=true)
+```
+
+The recipe says cook until the internal temperature is 165 degrees Fahrenheit. However, the food thermometer in the kitchen is calibrated in Celsius. Which value is closest?
+
+```{julia}
+#| echo: false
+choices = [68,70,72,74,76,78,80]
+v = 5/9*(165 - 32)
+_, answer = findmin(abs.(choices .- v))
+radioq(choices, answer; keep_order=true)
+```
+
+(Use $F = 9/5\cdot C + 32$ and $C = 5/9(F-32)$.)
+
+###### Question
+
 
 Which of the following is a valid `Julia` expression for
 
diff --git a/quarto/precalc/exp_log_functions.qmd b/quarto/precalc/exp_log_functions.qmd
index 1c583703..c8d416ef 100644
--- a/quarto/precalc/exp_log_functions.qmd
+++ b/quarto/precalc/exp_log_functions.qmd
@@ -8,7 +8,7 @@ This section uses the following add-on packages:
 
 ```{julia}
 using CalculusWithJulia
-using Plots
+using Plots; plotly()
 ```
 
 
diff --git a/quarto/precalc/functions.qmd b/quarto/precalc/functions.qmd
index ed7e9bfc..b0a3f3f0 100644
--- a/quarto/precalc/functions.qmd
+++ b/quarto/precalc/functions.qmd
@@ -8,6 +8,7 @@ This section will use the following add-on packages:
 
 ```{julia}
 using CalculusWithJulia, Plots
+plotly()
 ```
 
 
diff --git a/quarto/precalc/inversefunctions.qmd b/quarto/precalc/inversefunctions.qmd
index 36f2fb25..b650c3c5 100644
--- a/quarto/precalc/inversefunctions.qmd
+++ b/quarto/precalc/inversefunctions.qmd
@@ -9,6 +9,7 @@ In this section we will use these add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 ```
 
 ---
diff --git a/quarto/precalc/julia_overview.qmd b/quarto/precalc/julia_overview.qmd
index aa740b9d..ccb501ad 100644
--- a/quarto/precalc/julia_overview.qmd
+++ b/quarto/precalc/julia_overview.qmd
@@ -518,6 +518,7 @@ We must load `Plots` before we can plot (and it must be installed before we can
 
 ```{julia}
 using Plots
+plotly()  # optionally change the backend from the default
 ```
 
 With `Plots` loaded, we can plot a function by passing the function object by name to `plot`, specifying the range of `x` values to show, as follows:
diff --git a/quarto/precalc/logical_expressions.qmd b/quarto/precalc/logical_expressions.qmd
index cea00c9b..e0fec347 100644
--- a/quarto/precalc/logical_expressions.qmd
+++ b/quarto/precalc/logical_expressions.qmd
@@ -14,6 +14,7 @@ using CalculusWithJulia    # loads the `SpecialFunctions` package
 #| echo: false
 #| results: "hidden"
 using Plots
+plotly()
 nothing
 ```
 
diff --git a/quarto/precalc/plotting.qmd b/quarto/precalc/plotting.qmd
index 1cb42966..a574087e 100644
--- a/quarto/precalc/plotting.qmd
+++ b/quarto/precalc/plotting.qmd
@@ -9,6 +9,7 @@ This section will use the following packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 ```
 
 ```{julia}
@@ -74,8 +75,13 @@ using Plots
 
 :::
 
-The `plotly` backend is part of the `Plots` package, as is `gr`. Other backends require installation, such as `PyPlot` and `PlotlyJS`. We use `gr` in these notes, for the most part. (The `plotly` backend is also quite nice for interactive usage, but doesn't work as well with the static HTML pages.)
+Some backends require installation, such as `PyPlot` and `PlotlyJS`. We use `plotly` in these notes, for the most part, which is not the default, so requires an additional command to set the backend:
 
+```{julia}
+plotly()
+```
+
+(Certain graphics are produced with the `gr()` backend.)
 
 With `Plots` loaded, it is straightforward to graph a function.
 
diff --git a/quarto/precalc/polynomial.qmd b/quarto/precalc/polynomial.qmd
index 878db58a..02461b5d 100644
--- a/quarto/precalc/polynomial.qmd
+++ b/quarto/precalc/polynomial.qmd
@@ -9,6 +9,7 @@ In this section we use the following add-on packages:
 ```{julia}
 using SymPy
 using Plots
+plotly()
 ```
 
 ```{julia}
@@ -43,7 +44,7 @@ $$
 #| cache: true
 ##{{{ different_poly_graph }}}
 
-
+gr()
 anim = @animate for m in  2:2:10
     fn = x -> x^m
     plot(fn, -1, 1, size = fig_size, legend=false, title="graph of x^{$m}", xlims=(-1,1), ylims=(-.1,1))
@@ -53,6 +54,7 @@ imgfile = tempname() * ".gif"
 gif(anim, imgfile, fps = 1)
 caption = "Polynomials of varying even degrees over ``[-1,1]``."
 
+plotly()
 ImageFile(imgfile, caption)
 ```
 
@@ -87,7 +89,7 @@ $$
 #| echo: false
 #| cache: true
 ### {{{ lines_m_graph }}}
-
+gr()
 anim = @animate for m in  [-5, -2, -1, 1, 2, 5, 10, 20]
     fn = x -> m * x
     plot(fn, -1, 1, size = fig_size, legend=false, title="m = $m", xlims=(-1,1), ylims=(-20, 20))
@@ -96,7 +98,7 @@ end
 imgfile = tempname() * ".gif"
 gif(anim, imgfile, fps = 1)
 caption = "Graphs of y = mx for different values of m"
-
+plotly()
 ImageFile(imgfile, caption)
 ```
 
@@ -172,17 +174,17 @@ To illustrate, using the just defined `x`, here is how we can create the polynom
 
 
 ```{julia}
-𝒑 = -16x^2 + 100
+p = -16x^2 + 100
 ```
 
 That is, the expression is created just as you would create it within a function body. But here the result is still a symbolic object. We have assigned this expression to a variable `p`, and have not defined it as a function `p(x)`. Mentally keeping the distinction between symbolic expressions and functions is very important.
 
 
-The `typeof` function shows that `𝒑` is of a symbolic type (`Sym`):
+The `typeof` function shows that `p` is of a symbolic type (`Sym`):
 
 
 ```{julia}
-typeof(𝒑)
+typeof(p)
 ```
 
 We can mix and match symbolic objects. This command creates an arbitrary quadratic polynomial:
@@ -230,28 +232,28 @@ To illustrate, to do the task above for the polynomial $-16x^2 + 100$ we could h
 
 
 ```{julia}
-𝒑(x => (x-1)^2)
+p(x => (x-1)^2)
 ```
 
 This "call" notation takes pairs (designated by `a=>b`) where the left-hand side is the variable to substitute for, and the right-hand side the new value. The value to substitute can depend on the variable, as illustrated; be a different variable; or be a numeric value, such as $2$:
 
 
 ```{julia}
-𝒚 = 𝒑(x=>2)
+y = p(x=>2)
 ```
 
 The result will always be of a symbolic type, even if the answer is just a number:
 
 
 ```{julia}
-typeof(𝒚)
+typeof(y)
 ```
 
 If there is just one free variable in an expression, the pair notation can be dropped:
 
 
 ```{julia}
-𝒑(4) # substitutes x=>4
+p(4) # substitutes x=>4
 ```
 
 ##### Example
@@ -311,8 +313,8 @@ Evaluating a symbolic expression and returning a numeric value can be done by co
 
 
 ```{julia}
-𝐩 = 200 - 16x^2
-N(𝐩(2))
+p = 200 - 16x^2
+N(p(2))
 ```
 
 This approach is direct, but can be slow *if* many such evaluations were needed (such as with a plot). An alternative is to turn the symbolic expression into a `Julia` function and then evaluate that as usual.
@@ -323,7 +325,7 @@ The `lambdify` function turns a symbolic expression into a `Julia` function
 
 ```{julia}
 #| hold: true
-pp = lambdify(𝐩)
+pp = lambdify(p)
 pp(2)
 ```
 
@@ -356,7 +358,7 @@ Consider the graph of the polynomial `x^5 - x + 1`:
 plot(x^5 - x + 1, -3/2, 3/2)
 ```
 
-(Plotting symbolic expressions is similar to plotting a function, in that the expression is passed in as the first argument. The expression must have only one free variable, as above, or an error will occur.)
+(Plotting symbolic expressions with `Plots` is similar to plotting a function, in that the expression is passed in as the first argument. The expression must have only one free variable, as above, or an error will occur. This happens, as there is a `Plots` "recipe" for `SymPy` defined.)
 
 
 This graph illustrates the key features of polynomial graphs:
@@ -375,7 +377,7 @@ To investigate this last point, let's consider the case of the monomial $x^n$. W
 #| echo: false
 #| cache: true
 ### {{{ poly_growth_graph }}}
-
+gr()
 anim = @animate for m in  0:2:12
     fn = x -> x^m
     plot(fn, -1.2, 1.2, size = fig_size, legend=false, xlims=(-1.2, 1.2), ylims=(0, 1.2^12), title="x^{$m} over [-1.2, 1.2]")
@@ -384,7 +386,7 @@ end
 imgfile = tempname() * ".gif"
 gif(anim, imgfile, fps = 1)
 caption = L"Demonstration that $x^{10}$ grows faster than $x^8$, ... and $x^2$  grows faster than $x^0$ (which is constant)."
-
+plotly()
 ImageFile(imgfile, caption)
 ```
 
@@ -412,7 +414,7 @@ To see the roots and/or the peaks and valleys of a polynomial requires a judicio
 #| echo: false
 #| cache: true
 ### {{{ leading_term_graph }}}
-
+gr()
 anim = @animate for n in 1:6
     m = [1, .5, -1, -5, -20, -25]
     M = [2, 4,    5, 10, 25, 30]
@@ -427,7 +429,7 @@ end
 caption = "The previous graph is highlighted in red. Ultimately the leading term (\$x^4\$ here) dominates the graph."
 imgfile = tempname() * ".gif"
 gif(anim, imgfile, fps=1)
-
+plotly()
 ImageFile(imgfile, caption)
 ```
 
diff --git a/quarto/precalc/polynomial_roots.qmd b/quarto/precalc/polynomial_roots.qmd
index c99181d8..6507d21d 100644
--- a/quarto/precalc/polynomial_roots.qmd
+++ b/quarto/precalc/polynomial_roots.qmd
@@ -9,6 +9,7 @@ In this section we use the following add on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 using SymPy
 ```
 
diff --git a/quarto/precalc/polynomials_package.qmd b/quarto/precalc/polynomials_package.qmd
index 3cde9337..0d58b39c 100644
--- a/quarto/precalc/polynomials_package.qmd
+++ b/quarto/precalc/polynomials_package.qmd
@@ -9,6 +9,7 @@ This section will use the following add-on packages:
 ```{julia}
 import CalculusWithJulia
 using Plots
+plotly()
 using Polynomials
 using RealPolynomialRoots
 import SymPy # imported only: some functions, e.g. degree, need qualification
diff --git a/quarto/precalc/rational_functions.qmd b/quarto/precalc/rational_functions.qmd
index 92aacfc9..3244fb07 100644
--- a/quarto/precalc/rational_functions.qmd
+++ b/quarto/precalc/rational_functions.qmd
@@ -10,6 +10,7 @@ This section uses the following add-on packages:
 using CalculusWithJulia
 using SymPy
 using Plots
+plotly()
 import Polynomials
 using RealPolynomialRoots
 ```
diff --git a/quarto/precalc/transformations.qmd b/quarto/precalc/transformations.qmd
index 418caa4d..20b0e67a 100644
--- a/quarto/precalc/transformations.qmd
+++ b/quarto/precalc/transformations.qmd
@@ -9,6 +9,7 @@ In this section we will use these add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 ```
 
 ```{julia}
diff --git a/quarto/precalc/trig_functions.qmd b/quarto/precalc/trig_functions.qmd
index 752d8f65..3f5719a3 100644
--- a/quarto/precalc/trig_functions.qmd
+++ b/quarto/precalc/trig_functions.qmd
@@ -9,6 +9,7 @@ This section uses the following add-on packages:
 ```{julia}
 using CalculusWithJulia
 using Plots
+plotly()
 using SymPy
 ```
 
@@ -70,6 +71,7 @@ These definitions in terms of sides only apply for $0 \leq \theta \leq \pi/2$. M
 #| cache: true
 ## {{{radian_to_trig}}}
 
+gr()
 
 function plot_angle(m)
     r = m*pi
@@ -107,6 +109,7 @@ imgfile = tempname() * ".gif"
 gif(anim, imgfile, fps = 1)
 caption = "An angle in radian measure corresponds to a point on the unit circle, whose coordinates define the sine and cosine of the angle. That is ``(x,y) = (\\cos(\\theta), \\sin(\\theta))``."
 
+plotly()
 ImageFile(imgfile, caption)
 ```
 
@@ -293,7 +296,8 @@ The superposition of the two sine functions that `g2` represents could be viewed
 #| hold: true
 #| echo: false
 #| cache: true
-unzip(vs::Vector) = Tuple([[vs[i][j] for i in eachindex(vs)] for j in eachindex(vs[1])])
+gr()
+
 function makegraph(t, b₁,n₁, b₂=0, n₂=1)
 
     f₁ = x -> b₁*[sin(2pi*n₁*x), cos(2pi*n₁*x)]
@@ -350,7 +354,7 @@ imgfile = tempname() * ".gif"
 gif(anim, imgfile, fps = 5)
 
 caption = "Superposition of sines and cosines represented by an epicycle"
-
+plotly()
 ImageFile(imgfile, caption)
 ```
 
diff --git a/quarto/precalc/vectors.qmd b/quarto/precalc/vectors.qmd
index e7f53963..cd4e2c12 100644
--- a/quarto/precalc/vectors.qmd
+++ b/quarto/precalc/vectors.qmd
@@ -9,6 +9,7 @@
 #| results: "hidden"
 using CalculusWithJulia
 using Plots
+plotly()
 using Measures
 using LaTeXStrings
 
@@ -21,6 +22,8 @@ One of the first models learned in physics are the equations governing the laws
 ```{julia}
 #| hold: true
 #| echo: false
+
+gr()
 px = 0.26mm
 x0 = [0, 64]
 v0 = [20, 0]
@@ -45,7 +48,7 @@ function make_plot(t)
 
     xs, ys = [p[1] for p in xys], [p[2] for p in xys]
 
-    plt = Plots.plot(xs, ys, legend=false, size=fig_size, xlims=(0,45), ylims=(0,70))
+    plt = plot(xs, ys, legend=false, size=fig_size, xlims=(0,45), ylims=(0,70))
     plot!(plt, zero, extrema(xs)...)
     arrow!(xn(t), 10*unit(xn(t)), color="black")
     arrow!(xn(t), 10*unit(vn(t)), color="red")
@@ -74,6 +77,7 @@ end
 
 gif(anim, imgfile, fps = 1)
 
+plotly()
 ImageFile(imgfile, caption)
 ```
 
@@ -111,6 +115,8 @@ The magnitude of a vector comes from the distance formula applied to a line segm
 #| hold: true
 #| echo: false
 ## generic vector
+gr()
+
 p0 = [0,0]
 a1 = [4,1]
 b1 = [-2,2]
@@ -126,6 +132,7 @@ png(plt, imgfile)
 
 caption = "A vector and its unit vector. They share the same direction, but the unit vector has a standardized magnitude."
 
+plotly()
 ImageFile(imgfile, caption)
 ```
 
@@ -143,13 +150,13 @@ Two operations on vectors are fundamental.
 #| hold: true
 #| echo: false
 ## vector_addition_image
-
+gr()
 p0 = [0,0]
 a1 = [4,1]
 b1 = [-2,2]
 
 
-plt = Plots.plot(legend=false, size=fig_size)
+plt = plot(legend=false, size=fig_size)
 arrow!(p0, a1, color="blue")
 arrow!(p0+a1, b1, color="red")
 arrow!(p0, a1+b1, color="black")
@@ -160,6 +167,7 @@ png(plt, imgfile)
 
 caption = "The sum of two vectors can be visualized by placing the tail of one at the tip of the other"
 
+plotly()
 ImageFile(imgfile, caption)
 ```
 
@@ -167,7 +175,7 @@ ImageFile(imgfile, caption)
 #| hold: true
 #| echo: false
 ## vector_subtraction_image
-
+gr()
 p0 = [0,0]
 a1 = [-2,2]
 b1 = [4,1]
@@ -183,7 +191,7 @@ imgfile = tempname() * ".png"
 png(plt, imgfile)
 
 caption = "The difference of two vectors can be visualized by placing the tail of one at the tip of the other"
-
+plotly()
 ImageFile(imgfile, caption)
 ```
 
@@ -194,7 +202,7 @@ The concept of scalar multiplication and addition, allow the decomposition of ve
 #| hold: true
 #| echo: false
 ### {{{vector_decomp}}}
-
+gr()
 p0 = [0,0]
 aa = [1,2]
 bb = [2,1]
@@ -224,7 +232,7 @@ can be written uniquely as
 provided ``\vec{a}`` and ``\vec{b}`` are not parallel.
 
 """
-
+plotly()
 ImageFile(imgfile, caption)
 ```
 
@@ -248,6 +256,7 @@ If we are given $r$ and $\theta$, then the vector is $v = \langle r \cdot \cos(\
 #| hold: true
 #| echo: false
 ## vector_rtheta
+gr()
 p0 = [0,0]
 
 plt = plot(legend=false, size=fig_size)
@@ -266,6 +275,7 @@ A vector ``\langle x, y \rangle`` can be written as ``\langle r\cdot
 ``\theta``. The value ``r`` is a magnitude, the direction parameterized by
 ``\theta``."""
 
+plotly()
 ImageFile(imgfile, caption)
 ```