function [c,err,yc,p]=bisect(f,a,b,delta)
%Input - f is the function input as a string'f'
% - a and b are the left and right end points
% - delta is the tolerance
%0utput- c is the zero point
% - yc=f(c)
% - err is the error estimate for c
ya=feval(f,a);
yb=feval(f,b);
max1=1+round((log(b-a)-log(delta))/log(2));
p=0;
for k=1:max1
p=p+1;
c=(a+b) /2;
yc=feval(f,c);
if yc==0
a=c;
b=c;
elseif yb*yc>0
b=c;
yb=yc;
else
a=c;
ya=yc;
end
if b-a< delta
break
end
end
c=(a+b)/2;
err=abs(b-a);
yc=feval(f,c);
p;
end
