function p=steff(FUN, p0, TOL, MAX);
%%
%%   function p=p=steff(FUN, p0, TOL, MAX);
%%
%%   Steffensen"s mathod si to find a solution to 
%%   p=g(p) given an initial approximation p0;
%%
%%   
%%  INPUT:
%%  FUN  - the given function to be solved  which is entered as
%%           a string, enclosed in ' ', with the name the m-file of the
%%           function.
%%  p0 - the intiapproximation
%%  TOL - tolerance
%%  MAX - maximum number of interations
%%  OUTPUT:
%%  p- approximat solution p or massage of failure
i=1;
while i<=MAX
    i
    p1=feval(FUN, p0);
    p2=feval(FUN, p1);
    p=p0-(p1-p0).^2/(p2-2*p1+p0)
    if abs(p-p0)<TOL
       break; 
    end;
    i=i+1;
    p0=p;
end;
if i>MAX
     massage='mathod failed after N iteration'
end;
   