[p,Q]=steff(f,df,p0,delta,epsilon,max1)

function [p,Q]=steff(f,df,p0,delta,epsilon,max1)

%Input - f is the object function input as a string 'f'
%	    - df is the derivative of f input as a string 'df'
%      - p0 is the initial approximation to a zero of f
%      - delta is the tolerance for p0
%	    - epsilon is the tolerance for the function values y
%	    - max1 is the maximum number of iterations
%Output - p is the Steffensen approximation to the zero
%   	  - Q is the matrix containing the Steffensen sequence

% NUMERICAL METHODS: MATLAB Programs
%(c) 1999 by John H. Mathews and Kurtis D. Fink
%To accompany the textbook:
%NUMERICAL METHODS Using MATLAB,
%by John H. Mathews and Kurtis D. Fink
%ISBN 0-13-270042-5, (c) 1999
%PRENTICE HALL, INC.
%Upper Saddle River, NJ 07458

%Initialize the matrix R

R=zeros(max1,3);
R(1,1)=p0;

for k=1:max1
  for j=2:3

    %Denominator in Newton-Raphson method is calculated
    nrdenom=feval(df,R(k,j-1));

    %Conditional calculates Newton-Raphson approximations
    if nrdenom==0
      'division by zero in Newton-Raphson method'
       break		
    else
      R(k,j)=R(k,j-1)-feval(f,R(k,j-1))/nrdenom;
    end

    %Denominator in Aitkens Acceleration process is calculated
    aadenom=R(k,3)-2*R(k,2)+R(k,1);

    %Conditional calculates Aitkens Acceleration approximations
    if aadenom==0
      'division by zero in Aitkens Acceleration'
      break		
    else
      R(k+1,1)=R(k,1)-(R(k,2)-R(k,1))^2/aadenom;
    end
   end

   %Breaks out and ends program if division by zero occured
   if (nrdenom==0)|(aadenom==0)
     break
   end

   %Stopping criteria are evaluated; p and the matrix Q are determined	
    err=abs(R(k,1)-R(k+1,1));
    relerr=err/(abs(R(k+1,1))+delta);
    y=feval(f,R(k+1,1));
    if (err<delta)|(relerr<delta)|(y<epsilon)
       p=R(k+1,1);
       Q=R(1:k+1,:);
       break
    end
end