Example (E1) 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)-3*x(2)+x(1)^2*x(2)];

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

p =    0.15515*x1^2 + 0.056512*x1*x2 + 0.055778*x2^2;
zV = monomials(x,2:4);
Bis.flag = 0;
Bis.r1deg = 4;

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.923 	Beta UB = 0.923 
System 1: Num Stable = 9 	Num Unstable = 2 	Beta for Sims = 0.877 	Beta UB = 0.922 
System 1: Num Stable = 100 	Num Unstable = 2 	Beta for Sims = 0.877 	Beta UB = 0.922 
------------------End of simulations
------------------Begin search for feasible V
Try = 1 	 Beta for Vfeas = 0.877
------------------Found feasible V
Initial V (from the cvx outer bnd) gives Beta = 0.632
-------------------Iteration = 1 
Beta = 0.789 (Gamma = 0.770) 
-------------------Iteration = 2 
Beta = 0.805 (Gamma = 0.783) 
-------------------Iteration = 3 
Beta = 0.813 (Gamma = 0.790) 
-------------------Iteration = 4 
Beta = 0.818 (Gamma = 0.795) 
-------------------Iteration = 5 
Beta = 0.822 (Gamma = 0.799) 
-------------------Iteration = 6 
Beta = 0.826 (Gamma = 0.802) 
-------------------Iteration = 7 
Beta = 0.828 (Gamma = 0.804) 
-------------------Iteration = 8 
Beta = 0.831 (Gamma = 0.806) 
-------------------Iteration = 9 
Beta = 0.833 (Gamma = 0.808) 
-------------------Iteration = 10 
Beta = 0.835 (Gamma = 0.809) 
-------------------Iteration = 11 
Beta = 0.837 (Gamma = 0.811) 
-------------------Iteration = 12 
Beta = 0.838 (Gamma = 0.812) 
-------------------Iteration = 13 
Beta = 0.840 (Gamma = 0.813) 
-------------------Iteration = 14 
Beta = 0.841 (Gamma = 0.814) 
-------------------Iteration = 15 
Beta = 0.842 (Gamma = 0.814) 
-------------------Iteration = 16 
Beta = 0.842 (Gamma = 0.815) 
-------------------Iteration = 17 
Beta = 0.843 (Gamma = 0.816) 
-------------------Iteration = 18 
Beta = 0.844 (Gamma = 0.816) 
-------------------Iteration = 19 
Beta = 0.844 (Gamma = 0.817) 
-------------------Iteration = 20 
Beta = 0.845 (Gamma = 0.817) 

Extract the solution

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

    0.8447

betaUpper

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

    0.9224

Volume of the ROA estimate

Volume = vol/pi
Volume =

   14.9526