% Compare the following three methods:
% Newton's method,
% Secant method,
% Method of false position
TOL=1.e-8; N0=20;

p0 = 0.5; p1 = 0.5856;

if exist('method') == 0
    method = 'newton';
end

switch lower(method)
    
    case 'newton'
        
    i=1;
    while i <= N0 
        pause
        [y, yp] = ftosol(p0);
        p = p0 - y/yp;
        fprintf('approximation after %d iterations = %14.10f\n',i,p);
        if abs(p-p0) < TOL
            fprintf('found solution %14.10f\n',p);
            return;
        end
        i = i+1;
        p0=p;
    end
    fprintf('The method failed after %d iterations.\n',N0);
    
    case 'secant'
    
    i=2;
    q0 = ftosol(p0); q1 = ftosol(p1);
    while i <= N0
        pause
        p = p1 - q1*(p1-p0)/(q1-q0);
        fprintf('approximation after %d iterations = %14.10f\n',i,p);
        if abs(p-p1) < TOL
            fprintf('found solution %14.10f\n',p);
            return;
        end
        i = i+1;
        p0=p1; q0=q1; p1=p; q1=ftosol(p);
    end  
    fprintf('The method failed after %d iterations.\n',N0);
    
    case 'falsep'
        
    i=2;
    q0 = ftosol(p0); q1 = ftosol(p1);
    while i <= N0
        pause
        p = p1 - q1*(p1-p0)/(q1-q0);
        fprintf('approximation after %d iterations = %14.10f\n',i,p);
        if abs(p-p1) < TOL
            fprintf('found solution %14.10f\n',p);
            return;
        end
        i = i+1; q=ftosol(p);
        if q*q1 < 0
            p0=p1; q0=q1;
        end
        p1=p; q1=q;
    end  
    fprintf('The method failed after %d iterations.\n',N0);
    
    otherwise
    disp('method must be either Newton, Secant, or False Position.');

end
return