Example (E6) 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(1)+x(2)*x(3)^2;-x(2)+x(1)*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'*[0.2849    0.0356   -0.0008;
    0.0356    0.2888   -0.0014;
   -0.0008   -0.0014    0.2835]*x;
zV = monomials(x,2:4);
Bis.flag = 0;
Bis.r1deg = 4;

Generate the options

[roaconstr,opt,sys] = GetRoaOpts(f, x, zV, p, Bis);
opt.coordoptim.IterStopTol = 1e-4;
opt.display.roaest = 1;
opt.coordoptim.MaxIters = 30;
opt.sim.BetaInit = 2;
opt.sim.BetaMax = 2;

Solve the problem

outputs = wrapper(sys,[],roaconstr,opt);
------------------Beginning simulations
System 1: Num Stable = 100 	Num Unstable = 0 	Beta for Sims = 2.000 	Beta UB =  Inf 
------------------End of simulations
------------------Begin search for feasible V
Try = 1 	 Beta for Vfeas = 2.000
------------------Found feasible V
Initial V (from the cvx outer bnd) gives Beta = 1.401
-------------------Iteration = 1 
Beta = 1.612 (Gamma = 0.721) 
-------------------Iteration = 2 
Beta = 1.698 (Gamma = 0.763) 
-------------------Iteration = 3 
Beta = 1.739 (Gamma = 0.787) 
-------------------Iteration = 4 
Beta = 1.767 (Gamma = 0.804) 
-------------------Iteration = 5 
Beta = 1.789 (Gamma = 0.817) 
-------------------Iteration = 6 
Beta = 1.807 (Gamma = 0.827) 
-------------------Iteration = 7 
Beta = 1.821 (Gamma = 0.835) 
-------------------Iteration = 8 
Beta = 1.834 (Gamma = 0.844) 
-------------------Iteration = 9 
Beta = 1.843 (Gamma = 0.850) 
-------------------Iteration = 10 
Beta = 1.851 (Gamma = 0.855) 
-------------------Iteration = 11 
Beta = 1.858 (Gamma = 0.859) 
-------------------Iteration = 12 
Beta = 1.864 (Gamma = 0.863) 
-------------------Iteration = 13 
Beta = 1.869 (Gamma = 0.866) 
-------------------Iteration = 14 
Beta = 1.873 (Gamma = 0.868) 
-------------------Iteration = 15 
Beta = 1.877 (Gamma = 0.870) 
-------------------Iteration = 16 
Beta = 1.880 (Gamma = 0.873) 
-------------------Iteration = 17 
Beta = 1.884 (Gamma = 0.875) 
-------------------Iteration = 18 
Beta = 1.887 (Gamma = 0.877) 
-------------------Iteration = 19 
Beta = 1.890 (Gamma = 0.879) 
-------------------Iteration = 20 
Beta = 1.893 (Gamma = 0.880) 
-------------------Iteration = 21 
Beta = 1.896 (Gamma = 0.882) 
-------------------Iteration = 22 
Beta = 1.898 (Gamma = 0.883) 
-------------------Iteration = 23 
Beta = 1.900 (Gamma = 0.884) 
-------------------Iteration = 24 
Beta = 1.902 (Gamma = 0.886) 
-------------------Iteration = 25 
Beta = 1.904 (Gamma = 0.887) 
-------------------Iteration = 26 
Beta = 1.906 (Gamma = 0.888) 
-------------------Iteration = 27 
Beta = 1.908 (Gamma = 0.890) 
-------------------Iteration = 28 
Beta = 1.910 (Gamma = 0.891) 
-------------------Iteration = 29 
Beta = 1.912 (Gamma = 0.892) 
-------------------Iteration = 30 
Beta = 1.914 (Gamma = 0.893) 

Extract the solution

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

    1.9143

betaUpper

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

    20

Volume of the ROA estimate

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

   34.6967