Lyapunov theory region of attraction

V (y) = y Py

I have to find the region of attraction for this lyaponov function so I wrote this code but I couldn't finish it. (( I want to change the values of c to find the best ROA)

clear all

A = [-11 1 ; 3 -1];

Q = eye(2);

P = lyap(A',Q);

p11= P(1,1)

p12 = P(1,2)

p22 = P(2,2)

x=[-2:0.005:2];

y=x;

[X,Y]=meshgrid(x,y);

Z = -((X.^2)+(Y.^2))-((2*p11*X)+(2*p12*Y)).*(6+X).*(X.^2);

c = 0;

X1 = [];

X2 = [];

for i=1:1:size(y,2)

for j=1:1:size(x,2)

if Z(i,j)>=(c-0.01) && Z(i,j)<=(c+0.01);

X1=[X1;x(j)];

X2=[X2;y(i)];

end

end

end

f =figure

subplot(2,2,1)

plot(X1,X2)

clear X Y x1 x2 x y Z

c = 0.1;

X = [];

Y = [];

for x1= -2:0.005:2;

Delta = (2*p12*x1)^2-4*p22*((p11*x1^2)-c);

if Delta>0

x21 = (-2*p12*x1-sqrt(Delta))/(2*p22);

x22 = (-2*p12*x1+sqrt(Delta))/(2*p22);

X = [X;x1;x1];

Y = [Y;x21;x22];

end

end

subplot(2,2,2)

plot(X,Y)

c = 0.008;

Z = [];

W = [];

Answer 1

[Mo'ath Al-Hayek]

The final result shall be approximate to this