Bisection Method to Solve for a set of variables

I have two variables within a Bisection Method and I am trying to fill in a matrix that contains all the solutions to the equation. I can't seem to get the solutions of each unique set of variables to fill in the corresponding matrix. The two variables are M and T of length 24 and 13. so the solution matrix should look something like B = [ X1 of M1T1 X2 of M1T2 ....... [ X12 of M2T1 ......]

Is there something wrong with my loops?

for n = 1:length(Mach1)

    for m = 1:length(Theta1)

        xl = 0;

        xu = pi/6;

        tol = 0.001;

        f = inline('(tan(2*pi/180))-(2*(cot(x))*(((M)^2)*((sin(x))^2)-1))/(((M)^2)*(1.4+cos(2*x))+2)','x','M');

        for M = 2:25

            L = M - 1;

            for i = 2:1000

            xr=(xl+xu)/2;

                if f(xu,M)*f(xr,M)< 0

                   xl=xr;

                else

                   xu=xr;

                end

                if f(xl,M)*f(xr,M)<0

                    xu=xr;

                else

                    xl=xr;

                end

                xnew(1)=0;

                xnew(i)=xr;

                if abs((xnew(i)-xnew(i-1))/xnew(i))<tol

                    break

                end

                Beta1(L,1)=xr;

                xr=(xl+xu)/2;

            end