Guillaume, sorry that it has taken so long for someone to respond.
If your goal is to find beta and theta for each value of Ma1, you do that in the first two lines of the loop:
for i=1:length(Ma1)
beta=atand(1/Ma1(i)):1:90;
theta=atand(2*(cosd(beta)./sind(beta)).*(Ma1(i)^2*sind(beta).^2-1)./(Ma1(i)^2*(gamma+cosd(2*beta))+2));
plot(beta,theta)
hold on
end
This will plot a curve for each value of Ma1, if you provide a value for gamma. I added hold on so a curve is not wiped out when the next is plotted. I'm not sure what the rest of the code is for - maybe marking the maximum value of theta?