Next steps:
(0) (Relative Morse Theory) Geoffrey Mess' paper "Torelli groups of genus two and three surfaces" studies some relative Morse theory of the Abel-Jacobi period locus in the Siegel upper half spaces to deduce that the Torelli group (in genus two) is a free group on countably many generators. I thought his proof was very interesting, and tried to learn more, but hardly made progres...
(1) (Almost Complex Structures) if you're interested in symplectic topology, then Eliashberg-Cielebak's textbook "From Stein to Weinstein and back: Symplectic Geometry of Affine Complex Manifolds" has very interesting treatment of Morse theory, especially as related to almost-complex structures $J$ on symplectic manifolds $(M, \omega)$. I think this textbook eclipses Milnor's texts. Contains very elementary proof that "any $2n$-dimensional complex manifold has the homotopy type of an $n$-dimensional CW-complex". (Indeed the unstable manifold $W^+$ is totally lagrangian with respect to nondegenerate symplectic form $\omega=\omega_f$, and is therefore at most $n$-dimensional). Here $f$ is a real valued Morse function whose restriction to every $J$-invariant two-plane is subharmonic.
(2) Gradient flows to poles (where a potential function $f$ and its gradient $\nabla f$ diverges to $\pm \infty$) appears to have more applications to topology than the conventional gradient flow to zeros. Especially when attempting to strong deformation retract a noncompact source $X$ into a lower dimensional compact spine. Applying gradient flow to zeros requires a Lipschitz continuity-at-infinition condition on the deformation parameter. Here the Lowasiejiwicz inequality typically plays a decisive role in proving the continuity of the reparameterized gradient flow. The biggest problem with "gradient flow to zeros" is that the gradient flow slows down as it approaches its target. In my applications of optimal transport to algebraic topology, I find gradient flow to poles much more convenient, since the gradient enjoys a finite time blow up, and continuity of the reparameterized flow is immediate without any appeal to Lowasiejiwcz. Basically "gradient flow to zeros" is a soft landing, while "gradient flow to poles" accelerates into the target.
More specifically, I'm proposing that "gradient flow to poles" is important next step. And this occurs regularly in optimal transportation, as I describe next.
(3) (Optimal Transportation) Morse theory takes on new form in optimal transportation, where Morse theory plays a role in establishing the regularity/continuity and uniqueness of $c$-optimal transportation plans.
Consider a source probability space $(X, \sigma)$, target $(Y, \tau)$, and cost $c: X\times Y \to \mathbb{R}$. Kantorovich duality characterizes the $c$-optimal transport from $\sigma$ to $\tau$ via $c$-convex potential $\phi=\phi^{cc}$ on $X$ with $c$-transform $\psi=\phi^c$ on $Y$. Kantorovich says the $c$-optimal transport plan $\pi$ is supported on the graph of the $c$-subdifferential $\partial^c \phi$, or equivalently on the graph of $\partial^c \psi$.
The subdifferentials are characterized by the case of equality in $$-\phi(x)+\psi(y)\leq c(x,y).$$ Differentiating the case of equality with respect to $x$ and $y$ yields the equalities $$-\nabla_x \phi(x)=\nabla_x c(x,y)$$ and $$\nabla_y \psi(y)=\nabla_y c(x,y).$$ (R.J.McCann shows these equalities hold almost everywhere under general hypotheses on $c$). For example the (Twist) condition: If $Y\to T_x X$ defined by $y\mapsto \nabla_x c(x,y)$ is injective for every $x\in X$, then $$y=T(x):=\nabla_x c(x, \cdot)^{-1}(\nabla_x \phi(x))$$ defines a $c$-optimal Borel measurable map from $\sigma$ to $\tau:=T\#\sigma$.
Moreover the fibre $T^{-1}(y)$ can be characterized as the set of $x$ satisfying $\nabla_y\psi(y)=\nabla_y c(x,y)$ or $$\nabla_y [c(x,y)-\psi(y)]=0.$$ But observe that differentiating the $c$-Legendre Fenchel inequality a second time we are exlusively studying the global minimums of the potentials $y\mapsto c(x,y)-\psi(y)$, for every $x\in X$.
Using the usual Implicit Function theorem, the fibre $T^{-1}(y)$ is a smooth submanifold of $X$ if $D_x(\nabla_y c(x,y))$ is nondegenerate for every $x\in T^{-1}(y)$. If the target $(Y, \tau)$ is one-dimensional, this requires the function $x\mapsto \nabla_y c(x,y)$ to be critical point free for every $y\in Y$, and $x\in T^{-1}(y)$.
On most source manifolds $(X, \sigma)$ it is difficult to verify the nonexistence of critical points. If $X$ is compact and $c$ is continuous finite valued, then Morse theory (elementary calculus) forbids it. But we happily study costs $c$ with poles if the poles are the only critical values of $c$! For example, the (Twist) hypothesis can be rephrased as saying that the two pointed cross difference $$c_\Delta(x;y,y'):=c(x,y)-c(x,y')$$ is a critical point free function for all $y,y'$,$y\neq y'$ and $x$ on its domain. This cannot be satisfied on compact spaces unless poles are allowed.
(3.1) (Canonical Morse/Cost Functions?) We need distinguish generic and canonical. In my experience, I find generic functions very difficult to write down, or explore, or implement on Wolfram MATHEMATICA. Morse functions are known to be generic (in sense of Sard, Thom, etc.). But personally I prefer canonical Morse functions. Or from mass transport perspective, canonical costs $c$ whose derivatives $\nabla c$ are suitable Morse-type functions.
For example, if you want to study optimal transportation from a closed surface $\Sigma$ to the real line $Y=\mathbb{R}$ (or to circle or to graph), then one seeks an appropriate cost $c: \Sigma \times Y \to \mathbb{R}$ satisfying the above conditions, e.g. that $\frac{\partial c}{ \partial y}(x ,y)$ be critical point free in $x\in \Sigma$ for every $y\in \mathbb{R}$. This is forbidden by Morse theory if $\Sigma$ is compact and $c$ is everywhere finite. (In applications, we allow $c$ to have $+\infty$ poles. Then $\partial c/\partial y$ is possibly critical point free on its domain).
But what is a canonical cost $c: \Sigma \times \mathbb{R} \to \mathbb{R}$ which represents an interesting geometric transport from $\Sigma$ to $\mathbb{R}$? Here the source and target spaces $\Sigma$, $Y=\mathbb{R}$ have no interactions a priori, they are not even embedded in a common background space unless we suppose $Y\subset X$.