Example (E3) 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);-4*x(1)-3*x(2)-3*x(3)+x(1)^2*x(2)+x(1)^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.6352    0.2056    0.0745
    0.2056    0.2763    0.0778
    0.0745    0.0778    0.0402]*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-3;
opt.display.roaest = 1;
opt.sim.BetaMax = 1;
opt.sim.BetaInit = 1;

Solve the problem

outputs = wrapper(sys,[],roaconstr,opt);
------------------Beginning simulations
System 1: Num Stable = 17 	Num Unstable = 1 	Beta for Sims = 0.950 	Beta UB = 0.968 
System 1: Num Stable = 100 	Num Unstable = 1 	Beta for Sims = 0.950 	Beta UB = 0.968 
System 1: Num Stable = 128 	Num Unstable = 2 	Beta for Sims = 0.902 	Beta UB = 0.923 
------------------End of simulations
------------------Begin search for feasible V
Try = 1 	 Beta for Vfeas = 0.902
------------------Found feasible V
Initial V (from the cvx outer bnd) gives Beta = 0.447
-------------------Iteration = 1 
Beta = 0.568 (Gamma = 0.668) 
-------------------Iteration = 2 
Beta = 0.669 (Gamma = 0.760) 
-------------------Iteration = 3 
Beta = 0.712 (Gamma = 0.810) 
-------------------Iteration = 4 
Beta = 0.736 (Gamma = 0.842) 
-------------------Iteration = 5 
Beta = 0.754 (Gamma = 0.867) 
-------------------Iteration = 6 
Beta = 0.771 (Gamma = 0.891) 
-------------------Iteration = 7 
Beta = 0.781 (Gamma = 0.905) 
-------------------Iteration = 8 
Beta = 0.791 (Gamma = 0.919) 
-------------------Iteration = 9 
Beta = 0.796 (Gamma = 0.926) 
-------------------Iteration = 10 
Beta = 0.800 (Gamma = 0.932) 
-------------------Iteration = 11 
Beta = 0.804 (Gamma = 0.938) 
-------------------Iteration = 12 
Beta = 0.807 (Gamma = 0.943) 
-------------------Iteration = 13 
Beta = 0.810 (Gamma = 0.947) 
-------------------Iteration = 14 
Beta = 0.812 (Gamma = 0.951) 
-------------------Iteration = 15 
Beta = 0.815 (Gamma = 0.954) 
-------------------Iteration = 16 
Beta = 0.816 (Gamma = 0.956) 
-------------------Iteration = 17 
Beta = 0.818 (Gamma = 0.959) 
-------------------Iteration = 18 
Beta = 0.820 (Gamma = 0.962) 
-------------------Iteration = 19 
Beta = 0.821 (Gamma = 0.964) 
-------------------Iteration = 20 
Beta = 0.823 (Gamma = 0.966) 

Extract the solution

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

    0.8229

betaUpper

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

    0.9229

Volume of the ROA estimate

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

   21.4432