MATLAB: Find roots

asymptotefind all rootsfind all zerosfzeroperiodic function

Hi everyone.
I'm struggling with a pretty easy problem (I'm using Matlab since short!): I have to find the (1000 first) roots of the following function:
y = x*cot(x) + L -1
where c is a constant (in my case, usually between 1 and 100). I had to use a few tricks to avoid the undesired solutions that the fzero-function leads to. In the end, my code is the following (I called x beta):
function [beta, beta_null] = f_find_beta(L)
x0 = 0.001;
L = 10;
beta = zeros(1,1000);
beta_null = zeros(1,999);
%betaq = 0.1;
for q = 1:1:1000
      x = linspace( x0, x0 + pi );
      y = x .* cot(x) + L - 1;
      for p = 1:1:99
          if y( p+1 ) ./ y( p ) < 0
              z = fzero( @(x) x .* cot(x) + L - 1, x(p) );
              if abs( (z ./ pi) - (round(z ./ pi)) ) > 0.01 %&& (z - betaq) > 2.0
                  beta(q) = z;
                  %betaq = z;
              else end
          else end
      end
      x0 = x0 + pi ;
end
for n = 1:1:999
    beta_null(n) = beta(n+1) - beta(n);
end
end
beta_null is just a way for me to check my results more quickly. If you plot this vector as a function of its indices, you should get an almost straight line around pi.
It turns out that only the first 34-35 first roots can be calculated with this method without mistake. Then Matlab mistakes a "jump" of the function from +Inf to -Inf with a possible location of a root. I avoided the first mistakes with one of the if-loops (using a property of the x*cot(x)-function: its asymptotes are periodic (but not its roots) ) but the second trick did not improve my .m-file any further (see the "betaq" lines, with the % at the beginning since this method does not work).
Do you know how I could improve this file? In case this file seems to you to be a disaster, do you have another solution? (I also tried using the "findallzeros"-file that one can find really easily by typing "find roots matlab" or "find zeros matlab" in Google, but it does not work either).
Thanks for your help.
Tibo

Best Answer

The key to this problem is that there is one root in each interval [n*pi,(n+1)*pi]. I'm going to assume that L > 1, otherwise the region near the origin is a special case.
function beta = find_beta(L,nRoots)
f = @(x) x*cot(x)+L-1;
beta = zeros(nRoots,1);
% A small offset is needed to avoid the asymptotes.
smallOffset = 1e-12;
for i=1:nRoots
    beta(i) = fzero(f,[i*pi+smallOffset (i+1)*pi-smallOffset]);
end
Related Question