function [X,E,G]=fibonacci(f,a,b,tol,e)

%Input    - f, the object function 
%            - a, the left endpoint of the interval
%            - b, the right endpoint of the interval
%            - tol, length of uncertainty
%            - e, distinguishability constant
%Output - X, x and y coordinates of minimum
%Note this function calls the m-file fib.m

%If f is defined as an M-file function use the @ notation
% call X=fibonacci(@f,a,b,tol,e).
%If f is defined as an anonymous function use the
% call X=fibonacci(f,a,b,tol,e).

%Determine n
i=1;
F=1;
while F<=(b-a)/tol
    F=fib(i);
    i=i+1;
end

%Initialize values
n=i-1;
A=zeros(1,n-2);B=zeros(1,n-2);
A(1)=a;
B(1)=b;
c=A(1)+(fib(n-2)/fib(n))*(B(1)-A(1));
d=A(1)+(fib(n-1)/fib(n))*(B(1)-A(1));
k=1;
A(k)=a; B(k)=b; C(k)=c; D(k)=d;

%Compute Iterates
while k<=n-3
    fc=feval(f,c);
    fd=feval(f,d);
    if fc>fd
        A(k+1)=c;
        B(k+1)=B(k);
        c=d;
        d=A(k+1)+(fib(n-k-1)/fib(n-k))*(B(k+1)-A(k+1));
        C(k+1)=c; D(k+1)=d;
    elseif fc<fd
        A(k+1)=A(k);
        B(k+1)=d;
        d=c;
        c=A(k+1)+(fib(n-k-2)/fib(n-k))*(B(k+1)-A(k+1));
        C(k+1)=c; D(k+1)=d;
    else
        A(k+1)=c;
        B(k+1)=d;
        c=A(k+1)+(fib(n-k-2)/fib(n-k))*(B(k+1)-A(k+1));
        d=A(k+1)+(fib(n-k-1)/fib(n-k))*(B(k+1)-A(k+1));
        C(k+1)=c; D(k+1)=d;
    end
    k=k+1;
end

%Last iteration using distinguishability constant e
if feval(f,c)>feval(f,d)
    A(n-2)=c;
    B(n-2)=B(n-3);
    c=d;
    d=A(n-2)+(0.5+e)*(B(n-2)-A(n-2));
    C(n-2)=c; D(n-2)=d;
else
    A(n-2)=A(n-3);
    B(n-2)=d;
    d=c;
    c=A(n-2)+(0.5-e)*(B(n-2)-A(n-2));
    C(n-2)=c; D(n-2)=d;
end

%Output: Use midpoint of last interval for abscissa
if feval(f,c)>feval(f,d)
    a=c;b=B(n-2);ya=feval(f,a);yb=feval(f,b);
else
    a=A(n-2);b=d;ya=feval(f,a);yb=feval(f,b);
end
X=[(a+b)/2 feval(f,(a+b)/2)];
dp=abs(b-a);
dy=abs(yb-ya);
E=[dp dy];
G=[A' C' D' B'];