Note there is no need to use fzero. Roots is entirely adequate, and more accurate. Roots also recognizes there are TWO real solutions in general.
syms Area mach
>> gama = 1.4
>> Athroat = 650
(mach^2/6 + 5/6)^3/mach - Area/650
It is essentially just a polynomial. I've used syms just to make that clear.
Multiply by mach, fine as long as it is non-zero, and collect coefficients. Multiply by 216 too.
expand(((mach^2/6 + 5/6)^3/mach - Area/650)*mach)
mach^6 + 15*mach^4 + 75*mach^2 - (108*Area*mach)/325 + 125
So the polynomial coefficients as a function of Area are:
machcoeff = @(Area) [1 0 15 0 75 -108*Area*/325 125];
Now, we can use roots.
0.79428 + 3.8329i
0.79428 - 3.8329i
-1.9458 + 2.5675i
-1.9458 - 2.5675i
1.8863 + 0i
0.41673 + 0i
It finds two solutions that are not complex. The others are of no interest to us.
excludecomplex = @(vec) vec(imag(vec) == 0);
Now, put it all together:
Note there is no solution to the problem when Area is the same as Athroat.
0×1 empty double column vector
Just loop over the various values for Area.
Only you know which of those two solutions is meaningful, but as long as Area is greater than Athroat, two solutions seem to exist.
Areas = [660:20:1200];
Machs = NaN(2,length(Areas));
for n = 1:length(Areas)
M = excludecomplex(roots(machcoeff(Areas(n))));
Machs(:,n) = M;
I imagine one of those solutions is the one you are interested in seeing.
Of course, if you wanted to leave the other coefficients like Athroat and gama in there, you could have done that too.