Implicit system not working

Hello, I have troubles to run my code. I want to solve an implicit 2 variables-2 equations system, but it does not work. I have tried fzero fminsearch and fsolve but none of them works in my case. I think the problem comes from the way I've written my system but I cannot find any solutions and those already on the forum did not work here

Here is the main code

for i=1:length(Ma1)
    %1
    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));
    %2
    [thetam,I]=max(theta);
    B=[B thetam];
    betam=beta(I);
    T=[T betam];
    xmarkers=T;
    ymarkers =B;
    %3
    M1=Ma1(i)
    [B1,T1] = fzero(@myfun(beta,theta,M1),[60;0])
    %4
    plot(beta,theta,xmarkers,ymarkers,'b.')
    xlabel('\beta Angle')
    ylabel('\theta Angle')
    grid on
    axis([0 90 0 35])
end

The myfun function has been defined this way

function F = myfun(beta,theta,M1)
global gamma 
        F = [ ((gamma-1)*M1^2*sind(beta).^2+2)/(2*gamma*M1^2*sind(beta).^2+1-gamma)-sind(beta-theta).^2
              2*(cosd(beta)./sind(beta)).*(M1^2*sind(beta).^2-1)./(M1^2*(gamma+cosd(2*beta))+2)-tand(theta)];
end

I would like to find delta and theta for each value of Ma1 which is a 5 values vector If someone could help me with this, I would be very grateful Thanks

Answer 1

[Andrew Newell]

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?