function[p1,err,k,y]=newton(f1041,df1041,p0,delta,max1)
%  f是非線性函數。
%  df是f的微商。
%  p0是初始值。
%  delta是給定允許誤差。
%  max1是迭代的最大次數。
%  p1是牛頓法求得的方程的近似解。
%  err是p0的誤差估計。
%  k是迭代次數。
%  y=f(p1)
p0,feval('f1041',p0)
for k=1:max1
   p1=p0-feval('f1041',p0)/feval('df1041',p0);
       err=abs(p1-p0);
         p0=p1;
       p1,err,k, y=feval('f1041',p1)
         if(err<delta)| (y==0),
         break,
     end
     p1,err,k, y=feval('f1041',p1)
end