% defined functions
r = 1;
% r = 2;
g = @(x) x.*exp(r*(1-x));
f = @(x) x.*(1-exp(r*(1-x)));

% plot functions
m = 21;
x = 2*(0:(m-1))/(m-1);
fx = f(x);
gx = g(x);
plot(x,gx,'r',x,fx,'b',x,x,'k',x,zeros(size(x)),'c')

% Carry out bisection iteration
a = 0.1;
b = 1.5;
fa = f(a);
fb = f(b);
nmax = 1000000000;
ep = 1.0e-3;
de = 1.0e-3;
for n=1:nmax
    p = (a+b)/2;
    fp = f(p);
    if ((abs(fp) < ep) || (abs(b-a) < de*abs(p)))
        break;
    end
    if ((fa*fp) > 0)
        a = p; % sign(fa) = sign(fp)
    else
        b = p; % sign(fb) = sign(fp)
    end
end
if (n >= nmax)
    error('Bisection failed to converge after nmax iterations.');
end
disp(strcat('Bisektion p=',num2str(p),', iterations n=',num2str(n)));

% carry out fixed point iteration
p0 = 0.1;
for n=1:nmax
    p = g(p0);
    if ((abs(p-p0) < de*abs(p0)) | (abs(p-p0) < ep))
        break;
    end
    p0 = p;
end
if (n >= nmax)
    error('Fixed Point failed to converge after nmax iterations.');
end
disp(strcat('Fixed Point p=',num2str(p),', iterations n=',num2str(n)));

% results for r=1: (|g'(p)|=0)
% Bisektion p=0.99961, iterations n=9
% Fixed Point p=1, iterations n=6

% results for r=2: (|g'(p)|=1)
% Bisektion p=0.99961, iterations n=9
% Fixed Point p=0.9995, iterations n=2996981
