Is Lyapunov function the ultimate method to assess the stability analysis of a system?

The Lyapunov formalism is the standard approach to assessing the stability of classical perturbed nonlinear systems. As to why you don't see more emphasis placed upon it in physics education, I suspect it is because classical perturbation theory is just not covered as much in physics as it is in some areas of engineering.

This is my personal overly-simplified opinion: I would say that probably yes, Lyapunov function method is in most cases the best approach to stability in control theory, because it is a computationally friendly approach (and thus most often constructive method). Of course, there could be specific cases where there might be other methods.

The reason for this fact is that those working in control theory usually try to solve the converse problem of what those who work in physics try to solve.

Physics:

1. Starting from fundamental principles, build a model of the system dynamics in terms of differential equations.
2. Try to solve the differential equations in order to describe the time evolution of the system.
3. For a particular solution, describing the time evolution of a state of the system, figure out whether the solution is Laypunov (asymptotically) stable.

Control Theory:

1. Starting from fundamental principles, build a model of the system dynamics, without control, in terms of differential equations.
2. Given a desired time evolution of the system, find what additional terms (aka control) need to be added to the system's differential equations in order to force it to follow the desired time evolution.
3. Choose the added control so that the desired time evolution is a Lyapunov asymptotically stable solution of the controlled differential equations. This way you are guaranteeing that a large set of near-by initial configurations of the controlled system will eventually approach the desired behavior.

Lyapunov stability is a topological phenomenon and is complicated to justify in its raw form. Lyapunov functions are a constructive, computationally friendly and convenient method to test for Lyapunov stability and asymptotic stability. There is even a linerization method for checking Lyapunov stability of an equilibrium calculating the eigenvalues of a matrix. The difficult part is to find a Lyapunov function, which is often not a trivial endeavor.

Stability in physics: Start with equations of motion $to$ find a solution $to$ find a Lyapunov function (or another method) that guarantees stability of a solution (hard). Finding such a Lyapunov function is difficult. Frequently, a good choice is the Hamilton's function (Energy function) of some kind of near-by system (like in perturbation theory).

Stability in control theory: equations of motion without control plus a given desired time evolution $to$ pick a Lyapunov function (usually something simple, to your liking) $to$ find, using the Lyapunov function, control that makes the desired evolution a Lyapunov asymptotically stable solution of the controlled system.

So you see, in control theory, people choose the Lyapunov function and then look for the control. In physics, you cannot choose but you have to find a Lyapunov function, which is not easy.

(When I say physics, I also mean physics and mathematics.)

One of the advantages of the Lyapunov formalism, as opposed to other formalisms for analyzing stability, is the fact that it has the ability to draw global, rather than merely local, conclusions about the stability of the system. For example, one can compute the basin of attraction for a particular stable equilibrium using a properly-chosen Lyapunov function. As for the justification that it's "based on the energy of the system," this is a special case of a more general fact: you can use conservation laws to produce Lyapunov functions for physical systems. This is significant, because there is in general no systematic way to produce a Lyapunov function for an arbitrary ODE.