Example (E4) 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 x3;
x=[x1;x2;x3];
f = [-x(2);-x(3);-.915*x(1)+(1-0.915*x(1)^2)*x(2)-x(3)];

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 = x'*[12.5000   -8.1000    3.0000
   -8.1000   20.8000   -8.5000
    3.0000    -8.5000   13.4000]*x;
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 = 100 	Num Unstable = 0 	Beta for Sims = 5.000 	Beta UB =  Inf 
------------------End of simulations
------------------Begin search for feasible V
Try = 1 	 Beta for Vfeas = 5.000
------------------Found feasible V
Initial V (from the cvx outer bnd) gives Beta = 4.720
-------------------Iteration = 1 
Beta = 5.678 (Gamma = 1.052) 
-------------------Iteration = 2 
Beta = 6.024 (Gamma = 1.075) 
-------------------Iteration = 3 
Beta = 6.174 (Gamma = 1.084) 
-------------------Iteration = 4 
Beta = 6.248 (Gamma = 1.089) 
-------------------Iteration = 5 
Beta = 6.294 (Gamma = 1.092) 
-------------------Iteration = 6 
Beta = 6.318 (Gamma = 1.094) 
-------------------Iteration = 7 
Beta = 6.334 (Gamma = 1.096) 
-------------------Iteration = 8 
Beta = 6.345 (Gamma = 1.097) 
-------------------Iteration = 9 
Beta = 6.353 (Gamma = 1.098) 
-------------------Iteration = 10 
Beta = 6.359 (Gamma = 1.099) 
-------------------Iteration = 11 
Beta = 6.364 (Gamma = 1.099) 
-------------------Iteration = 12 
Beta = 6.367 (Gamma = 1.099) 
-------------------Iteration = 13 
Beta = 6.369 (Gamma = 1.099) 
-------------------Iteration = 14 
Beta = 6.371 (Gamma = 1.100) 
-------------------Iteration = 15 
Beta = 6.373 (Gamma = 1.100) 
-------------------Iteration = 16 
Beta = 6.375 (Gamma = 1.100) 
-------------------Iteration = 17 
Beta = 6.376 (Gamma = 1.100) 
-------------------Iteration = 18 
Beta = 6.377 (Gamma = 1.100) 
-------------------Iteration = 19 
Beta = 6.378 (Gamma = 1.100) 
-------------------Iteration = 20 
Beta = 6.379 (Gamma = 1.100) 

Extract the solution

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

    6.3790

betaUpper

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

    20

Volume of the ROA estimate

Volume = vol/(4*pi)*3
Volume =

    0.6636