dicho(f, ranges, epsi)

function [ranges, nbit] = dicho(f, ranges, epsi)
% [ranges, nbit] = dicho(f, ranges, epsi)
%      rsolution de f(x) = 0 par dichotomie.
%      Entres
%      f        nom de la fonction continue
%      ranges   matrices 2xn de rels dfinissant des intervalles
%               [ranges(1,i) ranges(2,i)] ou f n'a qu'une racine et tels que
%               f(ranges(1,i)).* f(ranges(2,i))< 0
%      epsi     erreur absolue desire sur les racines. (eps par dfaut)
%      Sorties
%      ranges   matrices 2xn de rels dfinissant des intervalles de
%               longueur au plus epsi contenant les racines et tels que
%               pour chaque i ranges(1,i) <= ranges(2,i).
%      nbit     retourne optionnellement le nombre d'itrations
%               effectu.
%
% La mthode est vectorise et calcule donc simultanment autant de
% racines de f(x)=0 que de colonnes dans ranges.
% la fonction intervs peut fournir une matrice de ranges.

if nargin <= 2, epsi = eps; end;

vals = feval(f, ranges);
if ~isempty(find(vals(1,:).*vals(2,:) > 0))
   error([ 'certains intervalles ne vrifient pas l''encadrement des racines']);
end;

ind  = find(vals(1,:)>= 0);
ranges([1,2], ind) = ranges([2, 1], ind); 

dx = ranges(2,:)-ranges(1,:);
x  = ranges(1,:);
fmids = feval(f, mean(x));
nbit = 0;

while 1
   ind = find( (abs(dx) > epsi) & (fmids ~= 0) );
   % ind contient les indices pour lesquels la convergence n'est pas atteinte
   if isempty(ind)   % toutes les racines ont t trouves
      break;
   end;
   nbit = nbit+1;
   dx(ind) = 0.5.*dx(ind);
   mids(ind) = x(ind) + dx(ind);
   fmids(ind) = feval(f, mids(ind));
   ind = find(fmids < 0);
   x(ind) = mids(ind);
end;

ranges = [x; x+dx];
ind  = find(ranges(1,:)>= ranges(2,:));  % remise en ordre
ranges([1,2], ind) = ranges([2, 1], ind);