function [k, p, err, P] = fixpt(g,p0,tol,max1)
% Input:   g is the iteration function input as a sring 'g' 
%          p0 is the initial guess for the fixed point
%          tol is the tolerance 
%          max1 is the maximum number of iterations
% Output:  k is the number of iterations that were carried out 
%          p is the approximation to the fixed point 
%          err is the error in the approximation
%          P contains the sequence {pn}
P(1)=p0;
for k=2:max1 
   P(k)=feval(g,P(k-1)); 
   err=abs(P(k)-P(k-1)); 
   relerr=err/(abs(P(k))+eps);
   p=P(k);
   if (err<: tol) | (relerr< tol), break;end
   end
if k==max1;
   disp('maximum number of iterations exceeded')
end
P=P';

% Fixed-Point Iteration Numerical Method for finding the x root of f(x) to make f(x) = 0

function [xR,err,n,xRV,errV,AFD1,AFD2] = FixedPointNM(AF,xi,ed)
% Inputs: with examples
% AF = anonymous function equation: AF = @(x) 1-((20^2)./(9.81*(((3*x)+((x.^2)/2)).^3))).*(3+x);
% xi = initial guess x = xR, where xR = x root: xi = 0.5;
% ed = desired approximate relative error = |(current - previous)/current|: ed = 0.01;
% Outputs
% xR = x root
% err = approximate relative error
% n = number of iterations
% xRV = x root vector
% errV = approximate relative error vector
% AFD1 = anonymous function 1st derivative
% AFD2 = anonymous function 2nd derivative
%% Computations
% Derivatives
AFSYM = sym(AF); % symbolic math function
AFD1SYM = diff(AFSYM); % symbolic math function 1st derivative
AFD2SYM = diff(AFD1SYM); % symbolic math function 2nd derivative
AFD1 = matlabFunction(AFD1SYM); % anonymous function 1st derivative
AFD2 = matlabFunction(AFD2SYM); % anonymous function 2nd derivative
% Fixed-Point Iteration Method
err = 1; % initial approximate relative error = 100%
k = 1; % initial counter
while ed < err % compares desired versus calculated error
    xR = xi+AF(xi); % x root using fixed point iteration method
    xRV(k+1) = xR; % stores the x root per iteration in a vector
    err = abs((xRV(k+1) - xRV(k))/xRV(k+1)); % approximate relative error
    errV(k+1) = err; % stores the error into a vector
    xi = xR; % new guess x = xR, where xR = x root
    k = k+1; % increase counter
end
n = k - 1; % number of iterations
xRV = xRV(2:end); % readjust x root vector
errV = errV(2:end); % readjust approximate relative error vector

function FixedPointIteration
%Fixed point iteration example.
%4x^2 = 10 - x^3 can be rewritten as x^2 = 0.25*(10 - x^3) 
%===> x = 0.5*(10- x^3)^0.5 if we take the postive root only.
%now we can use fixed point iteration.

i = 1;         
p0 = 1.5;           %our guess is p = 1.5.
error = 0.0000001;  %the precision we want to find our solution in.
N = 100;            %max number of iterations.

syms 'x' 
g(x) = 0.5*((10 - (x^3))^0.5); %i.e. what we are trying to solve

while i <= N
    p = g(p0);                  %the next p value
    if abs((p - p0)) < error    %stopping criterion using the required precision
        fprintf('Solution is %f \n', double(p))     %the output p is symbolic so need ot force it into a decimal
        return
    end
    
i = i + 1;   
p0 = p;         %update the value of i and p to continue the loop with

end
fprintf('Solution did not converge to within precision = %d in %d iterations \n', error, N)
end
 

% function f1=aitken(x,f,x1)
%  f1= aitken(x,f,x1)
%  f = corresponding function of x
%  f1= corresponding function of x1
   
    if ~isequal(length(x),length(f))
        disp('Error: Incomplete data points')
        f1=[];
        return
    end
    
    n=length(x); %length of available data
    D=x1-x;
    A(:,1)=f;
    for i=2:n
        for j=2:i
            A(i,j)=(D(j-1)*A(i,j-1)-D(i)*A(j-1,j-1))/(x(i)-x(j-1));
        end
    end
        
    f1=A(n,n);

function [p,Q] = steff(g,dg,p0,delta,epsilon,max1)
% Input:   g is the iteration function input as a sring 'g' 
%          dg is the derivative of g input as a string 'dg'
%          p0 is the initial guess for the fixed point
%          delta is the tolerance for p0
%          epsilon is the tolerance for the function values y
%          max1 is the maximum number of iterations
% Output:  Q is the matrix containing the Steffensen sequence 
%          p is the Steffensen approximation to the fixed point 

% Initialize the matrix R
R=zeros(max1,3);
R(1,1)=p0;
for k=1:max1 
    for j=2:3 
       % denominator in Newton method is calculated
       nrdenom=feval(dg,R(k,j-1));
       
       % calculate Newton approximations
       if nrdenom==0
          'division by zero in Newton method'
          break
       else
          R(k,j)=R(k,j-1)-feval(g,R(k,j-1))/nrdenom;
       end
    
    % denominator in Aitken's acceleration process 
       aadenom=R(k,3)-2*R(k,2)+R(k,1);
       
       % calculate Aitken's acceleration approximations
       if aadenom==0
          'division by zero in Aitken's acceleration'
          break
       else
          R(k+1,1)=R(k,1)-(R(k,2)-R(k,1))^2/aadenom;
       end
    
    end
    % end program if division by zero occured
    if (nrdenom==0)|(aadenom==0)  
       break
    end
      
    % stopping criteria are evaluated
    err=abs(R(k,1)-R(k+1,1));
    relerr=err/abs(R(k+1,1)+delta);
    y=feval(g,R(k+1,1)) ;
    if (err < delta)|(relerr < delta)|(y < epsilon)
       % p and the matrix Q are determined
       p=R(k+1,1);
       Q=R(1:k+1,:);
       break
    end
end

function  Complete

Complete