I completely agree with Luis's answer.   There is no general method to construct a Lyapunov function for ordinary differential equations.  However, for some specific system, one may use some method including the techniques from a view point of physical energy (the first integral in a term of math) to construct the Lyapunov function. Also most of the Lyapunov functions can only be used to deal with the local situation.  As for global cases or seeking the basin of attraction, the Lyapunov functions with nonnegative deviation could be taken into account.  

People above have pretty much described that there is no silver bullet for your question. Being said that, there are ways to find region of attraction (again numerically). The two topics that might be of interest to you would be "Sum of Squares" and "Contraction theory". "Sum of squares" is heavily used for finding region of attraction where as "Contraction theory" in one of its implementation can be used for the same. The basic idea is to numerically check Lyapunov function candidate and its Lie derivative around an equilibrium point.

Kindly refer to the link below to get some idea of what I am saying,

http://www-cdr.stanford.edu/dynamic/hcci_control/NOLCOS2004.pdf

The general answer to your question is given by Zubov's Theorem and related theorems. However, solving Zubov's theorem gives a partial differential equation which has to be solved for the Lyapunov function that defines the region of attraction exactly. This first order PDE has partial derivatives w.r.t all the states and hence is difficult to solve. This is still the simplest form to solve iteratively, starting from a known Lyapunov function from the control design. Most of the numerical methods to find Lyapunov functions are based on solving Zubov's pde. There are many mathematicians who have worked on this problem. Look up the work of Lars Grune, and Michael Malisoff and references therein to get a clear picture.

I agree with colleagues above. There is no general method to construct a Lyapunov function for ordinary differential equations. It depends on how the shape of your equations to see if it belongs to a system of literature by a change of variable for example. Another way to show the global stability of your steady state may be through analysis of nullclines. A good book for a introduction of Global Nonlinear Techniques is Differential equations, Dynamical Systems and  An Introduction to Chaos, Hirsch, Smale and Devaney.

Many thanks for all of your comments..

Could someone introduce a good reference material on Contraction theory and its application in determining the region of attraction in control engineering ?

The study of Nonlinear Dynamics has been centered on Fractals and Chaos Theory. Virgilio has given what I would say is a good reference because I recognize one of the authors; Robert L. Devaney. I have one of Devaney's books 'Chaos, Fractals and Dynamics', through which I began to understand what Chaos and Fractal Theories were about.

I have found that the book recommended by Virgilio is now available as a Third Edition. I am planning on buying a copy myself.

With regards to Contraction theory, I found by Googling:

https://books.google.com.au/books?id=L7Tw2yTKPZsC&pg=PA218&dq=Contraction+theory&hl=en&sa=X&ei=GGlxVceXBoSF8QWG7oGADg&redir_esc=y#v=onepage&q=Contraction%20theory&f=false

https://books.google.com.au/books?id=GzGo00c92tYC&pg=RA1-PA2&dq=Contraction+theory+control+systems&hl=en&sa=X&ei=v2pxVbb-PISymAXg54LwDQ&redir_esc=y#v=onepage&q=Contraction%20theory%20control%20systems&f=false

I'm not sure if you're familiar with Chaos and Fractal Theories. So I will give an introduction of my understanding of these branches of Nonlinear Dynamics.

A nonlinear system may Chaotically evolve over time, by taking one of two metastable trajectories. The trigger for which trajectory the nonlinear system ultimately follows can be a very small change in the meta-stability of the system. This is why a Chaotic system following one trajectory or the other, is referred to as existing in a metastable state because introducing another very small change will cause the trajectory to split or bifurcate into another two possible trajectories, and depending on the very small change, the system will adopt a new metastable trajectory defining the system's behaviour at that point.  Each trajectory may also evolve into two other separate trajectories. The pattern of these trajectories, may be referred to as Fractals.

The very small change to the system's metastable state defined by the current trajectory or Fractal parameter, may be thought of as perturbations imposed on the system.

Examples of Fractals are the bifurcation pattern of veins on the leaves of plants, or the pattern of blood vessels on an arm. Fractals do not have to be linear, they can be elliptical or circular, or a mixture such as coastlines. Other examples are the Mandelbrot set or the Julia set.

One interpretation of the term Contraction may be from Fractal mathematics, in which the smallest distance defining a Fractal as a multiple of the same basic pattern repeated in larger distances, is defined as contractive. Unless my understanding of contractive is flawed (from 'Mathematics of Digital Images; Creation, Compression, Restoration, Recognition' by S. G. Hoggar, 2006). Although the 'Dictionary of Mathematics' by Borowski and Borwein, 1989, defines contraction as 'n. 1. also called contraction mapping; a mapping on a METRIC SPACE that uniformly shrinks distances.'

Another concept from Chaos Theory, is the attractor. This is a region that some trajectories either orbit about, or are drawn into a central region or possibly points. The best example I can think of, is when making a cup of black coffee, and adding stevia or sugar, sometimes a frothy substance is formed and when stirring the cup, the froth closest to the boundary region between the coffee and the inner edge of the cup, will orbit about the center, and the froth closest to the center of the coffee whirlpool, will eventually make its way to the center. The froth, consisting of a large number of tiny bubbles, will result in the tiny bubbles tracing out the bifurcated trajectories of the liquid.

Another useful concept in Chaos Theory, is the limit cycle. A limit cycle is a type of attractor, that either exhibits an Attracting Limit Cycle, Repelling Limit Cycle or Saddle Cycle (one side has trajectory points repelled, the other side has trajectory points attracted to the attractor, or the other way around) within the Poincaré map function, for a two dimensional state space. The Attracting Cycle may be referred to as a stable limit cycle, and the Repelling Cycle may be referred to as an unstable limit cycle (from 'Chaos and Nonlinear Dynamics; An Introduction for Scientists and Engineers' by Robert C. Hilborn.

My assumption is that Contraction Theory may be a mathematical tool, allowing finding the smallest Fractal trajectory; although I'm not certain if this includes tools for finding trajectory points within the Stable Limit Cycle.

With regards to the Lyapunov function, in 'Introduction to Experimental Nonlinear Dynamics; A Case Study in Mechanical Vibration' by Lawrence N. Virgin, (2000), Lyapunov exponents are regarded as a key quantitative measure of Chaos. Trajectories can diverge rapidly, on the average and in a local sense this divergence takes place exponentially.

In the simplest Mathematical form: d(t)=d0*expλt

where d0 is the is the initial (small) separation of points on the attractor, d(t) is the distance between them at later times, and λ is the Lyapunov exponent (LE).

I hope the above information guides you to a better answer.

Excellent answers.

Still, there is no general solution for the general problem, so, along with reading, it is good to play a bit, or better more than just a bit, so you get some experience.

Before Lyapunov, dealing with stability of nonlinear systems simply meant not only to solve the nonlinear differential equation, but also to do it for all possible initial conditions, to make sure there are no diverging trajectories.

Lyapunov thought that, if a system is stable, then some energy-like function must increase if you go away from its assumed stable equilibrium point (or stable region) and decrease if you go towards it.

So, you may first try a simple Lyapunov function, such as V(x)=x^2 for a scalar system or V(x)=x1^2+x2^2+…+xn^2, where n is the number of variables. Now, if the derivative Vdot(x) of V(x) (and here the concept “along the trajectories” is the main idea) comes out to be negative definite, you are in business.

Along the trajectories is for the system xdot=f(x,t) and 

Vdot(x)=dv/dx dx/dt=dv/dx f(x,t). 

So, you substitute the particular f(x,t) and get the derivative  “along the trajectories.”

Because V was selected such that it cannot be negative, the derivative must end at Vdot=0, which is equivalent with x=0.

In more complex cases, Vdot is negative only outside some region and this gives you an idea about this region of (bounded) stability.

Even if the first simple V does not give you what you expected, the first attempts will give you more insight, so between the first V and the differential equation you may get closer to a proper V.

There are two excellent classical books on Lyapunov:

 [1] W. Hahn, Stability of Motion. New York, NY: Springer-Verlag, 1967.

[2] N. Rouche, P. Habets, and M. Laloy, Stability Theory by Lyapunov’s Direct Method. New York: Springer-Verlag, 1977.

Addenda:

I have another book 'Introduction to Nonlinear Differential and Integral Equations' by Harold T. Davis (reprint of the 1962 Edition). This book was published prior to the era of Chaos and Fractal theories, and in this book the author discusses the analysis of the roots of the integrated non-linear differential equations, and depending on the roots (real, imaginary or complex) and the signs of the roots, these define the phase trajectories of the differential equation under study, for seven different combinations of roots.

The book shows plots for Cases I to VII, and depending on the roots, the phase trajectories either converge towards a point, diverge from a point, or orbit about a point.

This book considers both nonlinear differential and partial differential equations; although the majority of the book is devoted towards the study of nonlinear differential equations.

Here is a link to the eBook version:

https://books.google.com.au/books?id=RgbWowrjKd4C&printsec=frontcover&dq=Introduction+to+Nonlinear+Differential+and+Integral+Equations&hl=en&sa=X&ei=RWFzVd_xJZfs8AW2yoGwAw&ved=0CB0Q6AEwAA#v=onepage&q=Introduction%20to%20Nonlinear%20Differential%20and%20Integral%20Equations&f=false

Whilst this book attempts to describe analytical techniques to find solutions to non-linear differential equations, these techniques are limited by the nature of Chaotic systems, to evolve radically when a very small perturbation is introduced, forcing the Chaotic system to bifurcate into a different metastable state. This is why the suggestion made by several respondents, to experiment with numerical solutions is so important.

As an interesting side note: Davis gives a practical application for the use of Limit Cycles, in finding the curve of best fit, for the nonlinear velocity graph for Cepheid Variable stars, from work done by W. S. Krogdahl, employing a generalized form of the Van Der Pol equation.

Cepheid Variable stars are stars which have regular pulsations. Astronomers use the brightness variations of CV stars in other galaxies, as "standard candles", to obtain reasonably accurate intergalactic distances up to 20 million light years. Cepheid Variable stars in our galaxy have had their distances verified by using parallax, up to the parallax limit of 65 light years from Earth's orbit, and with the launch of the Hipparcos satellite up to a parallax distance of 326 light years.

The Van Der Pol equation discussed in the book, was first used to describe the behaviour of the Triode oscillator in 1926.

HI!

You may take a look at the old Russian book by E. Barbashin entitled "Introduction to the theory of stability". Barbashin shows a way to constract Lyapunov functions for special classes of low-order systems. Even if it does not help you the book is worth reading in my opinion.

I found a supplier for E. Barbashin's book "Introduction to the theory of stability":

http://www.abebooks.com/9789001053451/Introduction-theory-stability-Barbashin-E-9001053459/plp

recommended by Dmitry.

I am planning on buying a copy for myself, to add to my own library, and to feed my curiosity to learn more about non linear dynamics.

As a side note, to numerically explore Chaos, Fractals, Lyapunov Exponents and Non linear dynamics, may I recommend using Mathcad, MATLAB, or their shareware counterparts SMath studio and SciLab or Octave?

For my Android based smart 'phone, I have installed a copy of a MATLAB like App called Addi, based on Octave:

https://play.google.com/store/apps/details?id=com.addi&hl=en

To generate plots, the AddiPlot App also needs to be downloaded and installed:

https://play.google.com/store/apps/details?id=com.addiPlot&hl=en

Re: Analytic Vs. Numerical solutions to Nonlinear Chaotic Fractal equations

Reading further Harold T. Davis' book "Introduction to Nonlinear Differential and Integral Equations", has given me cause to revise my previous comments on using numerical analyses approaches with a computer or programmable calculator; Davis writes in his Preface 'But the powers of the analytical method are not to be disregarded in favor of the powers of the machine. The solutions of nonlinear equations still possess singularities, which only the analytical method can discover and describe.'

The reason for Davis' comment, with which I agree, is because of the limits of numerical computation using Real numbers. Machines are limited in their representation of Real numbers because of the finite number of available bits used to represent these numbers. This leads to solutions when dealing with very small values, in which the digits farthest to the right of the decimal point hardly change at all, or for different parameters appear to converge to the same value; when intuitively we know that this should not be the case at all.

Case in point is when putting a number 2 in a calculator, then repeatedly pressing the square root key; eventually some calculators will give a value such as 1.00000000001 or some other calculators will show 1. Putting in a different number such as 3 will eventually give the same result. This is a limitation of using numerical approaches to solving either linear or non-linear equations; which should be considered when accepting the results of computerized numerical analyses.

The apparent convergence of different starting values to a similar point, or bounded region, is a major 'gotcha' of iterative or recursive numerical algorithms implemented in a computer with a sufficiently limited amount of bits or digits.

Perhaps using numerical approaches may be used as a guide to assist in finding analytical solutions, similar to how we use cars to get from one place to another quickly, and then we walk to get to a particular or specific place.

There is no  generally accepted approach. However some methods  exist, you can find them in D.R. Merkin, Introduction to the Theory of Stability, Springer-Verlag, 1997. N.G. Chetayev, The Stability of Motion, Pergamon Press, 1961.