The equation is a cubic in Unu^2, so the following uses roots to find solutions:
r = 3.88;
Vau = 43;
theta_u = 71.8;
gamma = 5/3;
Vsu = 47.45;
A = r*Vau^2*cos(theta_u)^2;
B = 2*r*Vsu^2/(r+1-gamma*(r-1));
C = r*Vau^2*sin(theta_u)^2;
D = (2*r - gamma*(r-1))/(r+1 - gamma*(r - 1));
% Let x = Unu^2
% (x - A)^2*(x - B) - C*x*(D*x - A) = 0
% (x/A - 1)^2*(x/A - B/A) - C/A*(x/A)*(D*x/A - 1) = 0
% Let y = x/A
% (y^2 - 2y + 1)*(y - B/A) - D*C/A*y^2 + C/A*y = 0
% Expand and collect like terms to get:
% y^3 -(2 + B/A + D*C/A)y^2 + (1 + 2*B/A + C/A)y - B/A = 0
p = [1; -(2 + B/A + D*C/A); (1 + 2*B/A + C/A); -B/A];
y = roots(p);
x = A*y; % reconstruct x
Unu = sqrt(x); % reconstruct Unu
disp(y)
disp(x)
disp(Unu)
% Check (f should be zero)
f = polyval(p,y);
disp(f)
