I suggest the following procedure only as approximation. Expand in Taylor series your differential equations at equilibrium point neighborhoods. The first order expansion is a linear approximation of your system in the equilibrium point neighborhoods and you can represent your system as,

dX/dt=Ax

the Lyapunov function can be,

V=XTPX

Where P is the solution of ATP+PA+Q=0  (with Q positive definite)

Analytically searching for Lyapunov functions is not a straight forward task. Because even if one function fails to prove stability, then one can not say that system is unstable. There might be another Lyapunov function for that. However numerically one can take the help of Sum Of Squares(SOS) programming for such issues.

Dear Ibrahim,

As mentionned already there is no general way of constructing Lyapunov functions for general systems of nonlinear differential equations. However, for systems with particular structures, such as linear or Hamiltonian (equiv. Lagrangian) for example, you can find certain generic strategies. I suggest you look at the books by Hassan Khalil, "Nonlinear Systems" and by F. Mazenc, M. Malissoff, "Constructions of Strict Lyapunov functions". Of course, given a nonlinear system, a first step would probably be to find a change of coordinates that can put your given system in a more convenient, simpler, form.

No way. There is no general method to find Lyapunov functions. It all depends on the structure of the particular system you are studying. If you would like to explore something a bit "algorithmic" you could perhaps pay a look to the integral averaging method. 

There is no general technique to find Lyapunov functions.  

What is the system that you will study?

The Lyapunov functions depends on the system.

If you are trying to carry out stability analysis of  a nonlinear system controlled by a fuzzy controller, an analytical way of determining a Lyapunov function whose partial derivative will be negative definite could be difficult if the number of fuzzy rules is large. The fewer the number of rules, the easier it will be to derive a Lyapunov function candidate. A good strategy is to use the Takagi-Sugeno modelling whose consequent parts are singleton. The search for bounds of stability will then be easier to compute and verify.