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) 2004 by John H. Mathews and Kurtis D. Fink%  Complementary Software to accompany the textbook:%  NUMERICAL METHODS: Using Matlab, Fourth Edition%  ISBN: 0-13-065248-2%  Prentice-Hall Pub. Inc.%  One Lake Street%  Upper Saddle River, NJ 07458%Initialize the matrix RR=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    endend