function [brackets,midpts] = nsection(fun,initinterval,n,tol)
% Find brackets for zeros of fun. 
% Generalizes bisection, looks for all bracket-intervals 
% Example:  
% fun=@sin;initint=[-10 10];n=11;
% [br,mid]=nsection(fun,initint,n);
%
if nargin==3
  tol=0.1;
end  
x=linspace(initinterval(1)-tol,initinterval(2)+tol,n);
y=fun(x);

sy=sign(y);     
d=diff(sy);        
sgnch=d~=0;        % Intervals with sign change
leftx=x(1:end-1);  % All left points
rightx=x(2:end);   % All right end points
leftbracket=leftx(sgnch);  
rightbracket=rightx(sgnch);
brackets=[leftbracket' rightbracket'];
midpts=sum(brackets,2)/2;
end

%{
clf
x=linspace(-10.1,10.1);y=sin(x);
plot(x,y);grid on;shg
hold on
lefts=br(:,1);rights=br(:,2);
plot(lefts,-0.1,'rx')
plot(rights,0.1,'bo')
stem(lefts,-0.1*ones(size(lefts)))
stem(rights,0.1*ones(size(rights)))
stem(mid,.2*ones(size(mid)))
%}