%Region of attraction example
%MCE/EEC 647/747
%Hanz Richter, PhD

%System:
%\dot{x}_1=-sin(x_1)cos(x_2) 
%\dot{x}_2=-cos(x_1)sin(x_2)

%Equilibrium points
%Unstable: (+/- k1*pi/2, +/- k2*pi/2) for k1,k2=0,1,2...
%Stable: (+/- k1*pi, +/- k2*pi) for k1,k2=0,1,2...

x10=pi;
x20=-pi;  %pick stable eq. point

%We use the quadratic Lyapunov function 
%V(x1,x2)=0.5((x1-pi)^2+(x2+pi)^2) which is p.s.d. near (pi,-pi)

%The derivative is 
%\dotV=-(x1-pi)sin(x1)cos(x2)-(x2+pi)cos(x1)sin(x2)

%First, since it's a 2D function,generate a contour plot for Vdot:

[X1,X2]=meshgrid([0:0.1:8],[-8:0.1:0]);
Vdot=-(X1-pi).*sin(X1).*cos(X2)-(X2+pi).*cos(X1).*sin(X2);
[c,h]=contour(X1,X2,Vdot);
clabel(c,h)
hold on
plot(pi,-pi,'+r')

%Now use a polar coordinate sweep to determine a radius of attraction

for r=0:0.1:5,
    for th=0:0.1:2*pi,
        x1=r*cos(th)+pi;
        x2=r*sin(th)-pi;  %absolute coords
        Vdot=-(x1-pi)*sin(x1)*cos(x2)-(x2+pi)*cos(x1)*sin(x2);
        if Vdot>0,
            r
            x1
            x2
            Vdot
            rectangle('Position',[pi-r -pi-r 2*r 2*r],'Curvature',[1,1]);
            axis equal
            error('Vdot>0. Terminating search')
        end
    end
end