function [p0,err,P] = fixpt(g,p0,tol,max1)
% [p0,err] = fixpt(g,p0,tol,max1)
% [p0,err,P] = fixpt(g,p0,tol,max1)
% Fixed point iteration.
% g is the function, input.
% p0 is the starting value, input.
% tol is the tolerance, input.
% max1 is the maximum number of iterations, input.
% p0 is the fixed point, output.
% err is the error estimate for p0, output.
% P is the vector of iterations, output.
P(1) = p0;
err = 1;
relerr = 1;
p1 = p0;
for k=1:max1,
p1 = feval(g,p0);
err = abs(p1-p0);
relerr = err/(abs(p1)+eps);
if (err<tol) | (relerr<tol), break; end
p0 = p1;
P(k+1) = p1;
end