%Bisection method
function [r,fr,logfc]=bisectwithreturn(fh,xa,xb)
if fh(xa)*fh(xb)>=0
    disp(['WARNING: BAD INITIAL GUESSES'])
    r=0;
    fr=0;
    logfc(1)=0;
    return
end
xc=(xa+xb)/2;
if fh(xc)==0
    disp(['Landed on a zero: stopping'])
    r=xc;
    fr=0;
    logfc(1)=0;
    return
end   
steps=0;
fcmin=1e-8;
fc=1;
while fc>fcmin && (steps<20)
     steps=steps+1;
     xold=xc;
     if fh(xc)*fh(xa)<0
         xb=xc;
     else
         xa=xc;
     end
     xc=(xa+xb)/2;
     fc=abs((xc-xold)/xc);
     logfc(steps)=log10(fc);
end %while loop
r=xc;
fr=fh(r);
end %funcition