Demonstrations used today:

• https://demonstrations.wolfram.com/FunctionsWithVerticalHorizontalAndObliqueAsymptotes/

WFR:

• https://resources.wolframcloud.com/FunctionRepository/resources/Asymptotes/
• https://resources.wolframcloud.com/FunctionRepository/resources/LucasCubic/

In the neat example section for the WFR LucasCubic you can find the non-trivial case used today for oblique asymptotes:

Demonstrations used Today:

• https://demonstrations.wolfram.com/PolynomialAndDerivative/
• https://demonstrations.wolfram.com/DerivativesALookAtGraphs/
• https://demonstrations.wolfram.com/CauchyMeanValueTheorem/
• https://demonstrations.wolfram.com/MeanValueTheorem/

Handy WFR function:

• https://resources.wolframcloud.com/FunctionRepository/resources/CurveAnalysis/

Here is an more involved example demonstrating that the parallel tangents can be found using numeric equation solving function:

f[x_, y_] := x^2 + x*Sin[4*y] + 3 y^2 - 7 + Cos[3*x];
pt = With[{x = RandomReal[{-2, 2}]}, {x, y /. FindRoot[f[x, y] == 0, {y, 1}]}];
slope = ImplicitD[f[x, y] == 0, y, x];
m = slope /. Thread[{x, y} -> pt];
pt2 = {x, y} /. FindRoot[{f[x, y] == 0, slope == m}, {x, -2}, {y, -1}];
m2 = slope /. Thread[{x, y} -> pt2];
ContourPlot[f[x,y]==0,{x,-4,4},{y,-4,4},Epilog->{InfiniteLine[pt,{1,m}],
    InfiniteLine[pt2,{1,m2}],
    {PointSize[0.02],Point[pt]},
    {PointSize[0.02],Point[pt2]}
}]

Having said that, here's a Lissajous code. Replace the first two lines in my original code with the following:

lj = {Sin[5 t + Pi/2], Sin[2 t]} // TrigToExp // ReplaceAll[t -> I*Log[z]];
f[x_, y_] = First@GroebnerBasis[{x, y} - lj, {x, y}, {z}]; 
boundsArray = {{-1, 1}, {-1, 1}};

You can replace {Sin[5 t + Pi/2], Sin[2 t]} with your favorite Lissajous parametrization — unless Mathematica chokes on it. In the above example, you see no tangents of slope 0. It's one of those problems I alluded to that I don't have time for; you have to put in checks and fixes for the various issues that might arise.

Thank you for the nice question! Here is the code for generating the plot shown on the book cover:

RemoveBackground[
 Rasterize[
  Plot3D[{Re[Sqrt[1/(Cos[5 x] Sin[5 y] - 1/2)]], 
    Im[Sqrt[1/(Cos[5 x] Sin[5 y] - 1/2)]], 
    Style[.01, White, Lighting -> {{"Ambient", White}}]}, {x, -1/3, 
    1.33}, {y, -1/3, 1.6}, Axes -> False, Boxed -> False, 
   BoundaryStyle -> None, Exclusions -> All, Mesh -> None, 
   PlotPoints -> 35, PlotRange -> {Automatic, Automatic, {0, 4}}, 
   ViewPoint -> {1.75, 2.25, 1.85}]]]

In order to understand the difference in the two plots, you could try setting Exclusions ->None for both of them.

The plot with PlotRange->{0, 4} is mostly smooth and setting Exclusions->All with a moderate setting of PlotPoints is enough to smooth out the singular behavior at the top.

The plot with PlotRange->{0, 8} is very irregular and using Exclusions->All with a much higher setting of PlotPoints (say PlotPoints->500) is required to get a reasonably smooth plot at the top, but it may still be difficult to get the same kind of smooth surface that is seen with PlotRange->{0,4}.