function [c,yc,err,P] = regula(f,a,b,delta,epsilon,max1)
% [c,yc,err] = regula(f,a,b,delta,epsilon,max1)
% [c,yc,err,P] = regula(f,a,b,delta,epsilon,max1)
% The Regula-Falsi method is used to locate a root.
% f is the function, input.
% a is the left endpoint, input.
% b is the right endpoint, input.
% delta is the tolerance for c, input.
% epsilon is the tolerance for yc, input.
% max1 is the maximum number of iterations, input.
% c is the root, output.
% yc is the function value, output.
% err is the error estimate for c, output.
% P is the is the matrix of endpoint iterations, output.
P = [a b];
ya = feval(f,a);
yb = feval(f,b);
if ya*yb > 0, break, end
for k=1:max1,
dx = yb*(b - a)/(yb - ya);
c = b - dx;
ac = c - a;
yc = feval(f,c);
if yc == 0,
break;
elseif yb*yc > 0,
b = c;
yb = yc;
else
a = c;
ya = yc;
end
P = [P;a b];
dx = min(abs(dx),ac);
if abs(dx) < delta, break, end
if abs(yc) < epsilon, break, end
end
err = abs(dx);