p0=1
delta=10^-5
epsilon=10^-5
for k=1:1000
   p1=p0-(3*f(p0))/(df(p0))
   err=abs(p1-p0);
   relerr=2*err/(abs(p1)+delta);
   p0=p1;
   y=f(p0)
   if (err<delta)|(relerr<delta)|(abs(y)<epsilon)
  break
  end
end

___________________________f.m
function y=f(x)
y=sin(x.^3);

__________________________df.m
function y=df(x)
y=(3*x.^2)*cos(x.^3);