function [c,err,yc]=bisect(f,a,b,delta) ya=feval(f,a); yb=feval(f,b); if ya*yb>0,break,end max1=-1+round((log(b-a)-log(delta))/log(2)); for k=1:max1 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; y