% regfalsi.m: Does false position method on given function.
% Usage: regfalsi('f', eps, a, b).  eps is tolerance, [a,b] bracket.
function c = regfalsi(f, eps, a, b)

fa = f(a); fb = f(b);

if fa*fb > 0
    disp('[a,b] doesn''t bracket, moron');
    c = NaN;
    return
end

k = 1;

while b-a>eps && k < 100000 %Dumb termination criterion, but OK for now.
    
    c = (a*fb-b*fa)/(fb-fa);
    fc = f(c);
    disp([k,a,c,b]);
    
    if fc*fa < 0 %[a,c] is a bracket
        b = c; fb = fc;
    else %[c,b] is the bracket
        a = c; fa = fc;
    end
    
    k = k+1;
    
end