ROA of SMIB system

[Vurgit]

We have tested this on a SMIB system, but something seems to go wrong. One with polynomialization returns ROA with volume = 0 and one with taylor expansion returns volume = 1.7 but both of them cannot plot the ROA. We are wondering why this will happen and codes are attached. Thank you for your patience and time if you could help us with this.
`clear
syms w theta
syms z1 z2
H = 0.5;
D = 0.1;
X = 0.1;
E = 1;
V = 1;
Pm = 0.9;
f = [
w;
(Pm - E * V * sin(theta) / X - D * w) / 2 / H;
];
x = [w; theta];
x0 = [0; asin(X * Pm)];
% sinθ = z1, cosθ = z2
f_tayl = taylor(f, x, x0, 'Order', 4);
% x_recast = [w; z1; z2];
% x0_recast = [0; sin(x0(2)); cos(x0(2))];

eq = x0;
dev = [1; 1];
% ang_ind = [2, 3];
SMIB = SOStab(eq, dev);
sysDyn = sym2SOStab('SMIB', f_tayl, x);
SMIB.dynamics = eval(sysDyn);
[vol, vc, wc] = SMIB.SoS_out(4,5,0.1);
vopow = value(SMIB.vcoef_outer);
wopow = value(SMIB.wcoef_outer);
SMIB.plot_roa(1, 2, 'outer');`

`clear
syms w theta
syms z1 z2
H = 0.5;
D = 0.1;
X = 0.1;
E = 1;
V = 1;
Pm = 0.9;
f = [
w;
(Pm - E * V * sin(theta) / X - D * w) / 2 / H;
];
x = [w; theta];
x0 = [0; asin(X * Pm)];
% sinθ = z1, cosθ = z2
f_recast = [
w * z2;
w * (-z1);
(Pm - E * V * z1 - D * w) / 2 / H;
];
x_recast = [w; z1; z2];
x0_recast = [0; sin(x0(2)); cos(x0(2))];

eq = x0_recast;
dev = [1; 1; 1;];
ang_ind = [2, 3];
SMIB = SOStab(eq, dev, ang_ind);
sysDyn = sym2SOStab('SMIB', f_recast, x_recast);
SMIB.dynamics = eval(sysDyn);
[vol, vc, wc] = SMIB.SoS_out(10,5,0.1);
vopow = value(SMIB.vcoef_outer);
wopow = value(SMIB.wcoef_outer);
SMIB.plot_roa(1, [2, 3], 'outer');`

function sysDyn = sym2SOStab(sysName, f, x) if ~isa(sysName, 'string') sysName = string(sysName); end f_str = string(f); x_str = string(x); [a, b] = size(x); nx = numel(x); dots = repmat(".", a, b); sysPre = join([repmat(sysName, a, b), dots], ""); % xPre = "x"; xNew_str = join([repmat("x(", a, b), join([string((1 : nx).'), repmat(")", a, b)]) ], ""); f_str_new = replace(f_str, x_str, join([sysPre, xNew_str], "")); sysDyn = f_str_new; sysDyn = join([sysDyn, repmat(";", nx, 1)]); sysDyn = join(sysDyn.'); sysDyn = join(["[", sysDyn, "]"]); end