p0=1
delta=10^(-5)
epsilon=10^(-5)
for k=1:1000
   p1=p0-(f(p0)*df(p0))/(df(p0)*df(p0)-f(p0)*ddf(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


__________________________ddf.m
function y=ddf(x)
y=6*x.*cos(x.^3)-(9*x.^4)*sin(x.^3);