Roots of Equations - Numerical Methods in Engineering with MATLAB - page 152

if f1*f2 > 0.0
error('Root is not bracketed in (a,b)')
end
x3=0.5*(a+b);
% Beginning of iterative loop.
fori=1:30
f3 = feval(func,x3);
if abs(f3) < tol
root = x3; return
end
% Tighten brackets (a,b) on the root.
iff1*f3<0.0;b=x3;
else;a=x3;
end
if (b - a) < tol*max(abs(b),1.0)
root = 0.5*(a + b); return
end
% Try quadratic interpolation.
denom = (f2 - f1)*(f3 - f1)*(f2 - f3);
numer = x3*(f1 - f2)*(f2 - f3 + f1)...
+ f2*x1*(f2 - f3) + f1*x2*(f3 - f1);
% If division by zero, push x out of bracket
% to force bisection.
ifdenom==0;dx=b-a;
else; dx = f3*numer/denom;
end
x=x3+dx;
% If interpolation goes out of bracket, use bisection.
if(b-x)*(x-a)<0.0
dx=0.5*(b-a);x=a+dx;
end
%Letx3<--x&choosenewx1,x2sothatx1<x3<x2.
ifx<x3
x2=x3;f2=f3;
else
x1=x3;f1=f3;
end
x3=x;
end
root = NaN;