Need a better method to plot this function

Essentially I am trying to find a better way of solving this:

    clear all; close all; home
N=10;
M=100;
L=400;
ki=zeros(L-1,floor(2/3*M+1));
kr=zeros(L-1,M-floor(2/3*M+1)+1);
options=optimset('MaxIter', 1e8, 'TolFun', 1e-6);
for i=1:floor(2/3*M+1)
    p=pi/M*(i-1);
    g_p=2*cos(p/2);
    for j=1:L-1
        g=@(x) sin(x*N)+g_p*sin(x*(N+1)); 
        [ki(j,i), feval1, exitflag1(j,i)]=fsolve(g,pi/L*j,options);
        if (exitflag1(j,i)~=1)
            ki(j,i)=0;
        end
    end
end
for i=(floor(2/3*M+1)+1):(M+1)
    m=i-floor(2/3*M+1);
    p=pi/M*(i-1);
    g_p=2*cos(p/2);
    for j=1:L-1
        g=@(x) sin(x*N)+g_p*sin(x*(N+1)); 
        [kr(j,m), feval2, exitflag2(j,m)]=fsolve(g,pi/L*j,options);    
        if (exitflag2(j,m)~=1)
            kr(j,m)=0;
        end
   end
end

which is nothing more than a set of nonlinear equations. The problem is that it takes too long for those given parameters. I know the system has the following properties: if i<=2/3*M+1 then we have N different real solutions; if i>2/3*M+1 then we have N-1 different real solutions. In both cases, the solutions have to lie between 0 and pi, without the boundary as possible solutions. The next piece of code just gets rid of the invalid solutions and sorts them in crescent order.

 for j=1:floor(2/3*M+1)
    ki(:,j)=sort(ki(:,j));  
 end
for j=(floor(2/3*M+1)+1):(M+1)
    m=j-floor(2/3*M+1);
    kr(:,m)=sort(kr(:,m));
end   
Ki=zeros(L,floor(2/3*M+1));
Kr=zeros(L-1,M-floor(2/3*M+1)+1);
for j=1:floor(2/3*M+1)
    m=0;
    for i=1:L-1
        if ((abs(Ki(:,j)-ki(i,j))>1e-5)&(ki(i,j)>0)&(ki(i,j)<3.1416))
        m=m+1;
        Ki(m,j)=ki(i,j);
        end
   end
end
for j=(floor(2/3*M+1)+1):(M+1)
    n=j-floor(2/3*M+1);
    m=0;
    for i=1:L-1
        if ((abs(Kr(:,n)-kr(i,n))>1e-5)&(kr(i,n)>0)&(kr(i,n)<3.1416))
        m=m+1;
        Kr(m,n)=kr(i,n);
        end
   end
end

And effectively, we get N different values for each column of Ki and N-1 different values for each column of Kr. I am seeking a faster and neater way of solving this problem making sure that for a given N (the actual problem requires N=30), it gets all possible solutions. Thanks.

Note: More information about the system I am dealing with can be found here http://iopscience.iop.org/article/10.1088/1468-6996/11/5/054504/meta, appendix B, look for the dispersion relation.

Answer 1

[Alan Weiss]

Without understanding what you wrote in any detail, it seems to me that you are looking at a large number of real roots of a real function. If that is true, perhaps you should try chebfun.

Alan Weiss

MATLAB mathematical toolbox documentation

[JLPiqueres]

Without understanding what I wrote you showed me a nice shortcut. Thanks!