Example (E2) from the paper (*)

(*) Topcu, Packard, Seiler, "Local stability analysis using simulations and sum-of-squares programming," Automatica, 2008.

Form the vector field

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

p, the basis vector for V and the iteration options (if not equal to the default values) need to be specified before generating the options structures

zV = monomials(x,2:4);
Bis.flag = 1;
Bis.s2deg = 2;
p = x'*[0.9937    0.2283
    0.2283    0.8907]*x;

Generate the options

[roaconstr,opt,sys] = GetRoaOpts(f, x, zV, p, Bis);

Modify the options

opt.coordoptim.IterStopTol = 1e-4;
opt.display.roaest = 1;

Solve the problem

outputs = wrapper(sys,[],roaconstr,opt);
------------------Beginning simulations
System 1: Num Stable = 0 	Num Unstable = 1 	Beta for Sims = 0.840 	Beta UB = 0.840 
System 1: Num Stable = 100 	Num Unstable = 1 	Beta for Sims = 0.840 	Beta UB = 0.840 
------------------End of simulations
------------------Begin search for feasible V
Try = 1 	 Beta for Vfeas = 0.840
------------------Found feasible V
Initial V (from the cvx outer bnd) gives Beta = 0.547
-------------------Iteration = 1 
Beta = 0.735 (Gamma = 0.796) 
-------------------Iteration = 2 
Beta = 0.787 (Gamma = 0.830) 
-------------------Iteration = 3 
Beta = 0.787 (Gamma = 0.830) 

Extract the solution

[V,beta,gamma,p,multip,betaUpper] = extractSol(outputs);
beta
beta =

    0.7871

betaUpper

[vol,volstd] = pvolume(V,gamma,[],10000);
betaUpper =

    0.8396

Volume of the ROA estimate

Volume = vol/pi
Volume =

    0.9997