pvar x1 x2;
x = [x1;x2];
x1dot = -x2;
x2dot = x1+(x1^2-1)*x2;
f = [x1dot; x2dot];
% Lyap func from linearization
Q = eye(2);
Vlin = linstab(f,x,Q);
% maximize gamma
% subject to:
% {Vlin<=gamma} in {Vdot<0} U {x=0}
z = monomials(x, 1:2 );
L2 = 1e-6*(x'*x);
Vdot = jacobian(Vlin,x)*f;
[gbnds,s] = pcontain(Vdot+L2,Vlin,z);
Gamma = gbnds(1);

% Compute and plot the limit cycle
[xtraj,xconv] = psim(f,x,[2;-2],-20);
[xtraj,xconv] = psim(f,x,xtraj{2}(end,:)',-20);
plot(xtraj{2}(:,1),xtraj{2}(:,2),'b','linewidth',2);
hold on;

% Plot the computed sublevel set
[C,h] = pcontour(Vlin, Gamma,[-3 3 -3 3]);
grid on;
set(h,'linewidth',2,'color','g');

% Plot the sublevel set of Vdot
Vdot = jacobian(Vlin,x)*f+L2;
[C,h] = pcontour(Vdot, 0 ,[-3 3 -3 3]);
set(h,'linewidth',2,'color','r');

% Sample the set {Vlin=Gamma} and simulate the systems from those points
[xin,xon] = psample(Vlin-Gamma,x,zeros(2,1),25);
[xtraj,xconv] = psim(f,x,xon,10);
for i1 = 1:size(xon,2)
    plot(xtraj{i1,2}(:,1),xtraj{i1,2}(:,2),'r');
    plot(xtraj{i1,2}(1,1),xtraj{i1,2}(1,2),'k*');
end

% Plot a few more sublevel sets of V
[C,h] = pcontour(Vlin, [0.2:0.2:0.8]*Gamma,[-3 3 -3 3]);
grid on;
set(h,'linewidth',2,'color','b');


% Pick a point in {Vdot < 0} but not in {V<=Gamma}
xd = [2.5;0];
[xtraj,xconv] = psim(f,x,xd,1);
plot(xtraj{2}(:,1),xtraj{2}(:,2),'m','linewidth',2);
hold on;