pvar x1 x2;
x = [x1;x2];
x1dot = -x2;
x2dot = x1 + (x1^2-1)*x2;
f = [x1dot; x2dot];

% 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),'k','linewidth',2);
hold on;

% Create shape function and monomials vectors
p = x'*x;
zV = monomials( x, 2:4 );
z1 = monomials( x, 0:1 );
z2 = monomials( x, 1:2 );
L2 = 1e-6*(x'*x);
% Initialize Lyapunov Function
V = linstab(f,x);
% Run V-s iteration
opts.L2 = L2;
for i1=1:25;
    % gamma step
    Vdot = jacobian(V,x)*f;
    [gbnds,s2] = pcontain(Vdot+L2,V,z2,opts);
    gamma = gbnds(2);
    co = 'b';
    if i1 == 1 co = 'g'; end
    if i1 == 25 co = 'r'; end
    [Ca,ha] = pcontour(V,gamma,[-3 3 -3 3]);
    set(ha,'color',co);
    hold on;
    % beta step
    [bbnds,s1] = pcontain(V-gamma,p,z1,opts);
    beta = bbnds(1)
    % V step (then scale to roughly normalize)
    if i1~=25
        V = roavstep(f,p,x,zV,beta,gamma,s1,s2,opts);
        V = V/gamma;
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%
% repeat with 6th degree V
close all;
[Ca,ha] = pcontour(V,gamma,[-3 3 -3 3]);
set(ha,'color','m','linewidth',2);
hold on;
clear all;

pvar x1 x2;
x = [x1;x2];
x1dot = -x2;
x2dot = x1 + (x1^2-1)*x2;
f = [x1dot; x2dot];

% 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),'k','linewidth',2);
hold on;

% Create shape function and monomials vectors
p = x'*x;
zV = monomials( x, 2:6 );
z1 = monomials( x, 0:2 );
z2 = monomials( x, 1:2 );
L2 = 1e-6*(x'*x);
% Initialize Lyapunov Function
V = linstab(f,x);
% Run V-s iteration
opts.L2 = L2;
for i1=1:25;
    % gamma step
    Vdot = jacobian(V,x)*f;
    [gbnds,s2] = pcontain(Vdot+L2,V,z2,opts);
    gamma = gbnds(2);
    co = 'b';
    if i1 == 1 co = 'g'; end
    if i1 == 25 co = 'r'; end
    [Ca,ha] = pcontour(V,gamma,[-3 3 -3 3]);
    set(ha,'color',co);
    hold on;
    % beta step
    [bbnds,s1] = pcontain(V-gamma,p,z1,opts);
    beta = bbnds(1)
    % V step (then scale to roughly normalize)
    if i1~=25
        V = roavstep(f,p,x,zV,beta,gamma,s1,s2,opts);
        V = V/gamma;
    end
end