編寫Newton迭代法的M文件NewtonIterate.如下：
function s=NewtonIterate(x,eps)
%Newton迭代法求非線性方程組
%x為迭代初值，eps為允許誤差值
if nargin==1
    eps=1.0e-6;
elseif nargin<1
    error
    return
end  
x1=fx1(x);%非線性方程組
x2=-dfx1(x);%非線性方程組導數
x3=inv(x2);
x0=x3*x1';
while norm(x0)>=eps %循環迭代
      x=x0'+x;
      x1=fx1(x);
      x2=-dfx1(x);
      x3=inv(x2);
      x0=x3*x1';
end
s=x0'+x;
return
例： x1.^2-10x1+x2.^2+8=0
     x1*x2.^2+x1-10x2+8=0
用MATLAB實現，編寫非線性方程組的M文件fx1.m如下所示：
   function y=fx1(x)
   y(1)=x(1)*x(1)-10*x(1)+x(2)*x(2)+8;
   y(2)=x(1)*x(2)*x(2)+x(1)-10*x(2)+8;
   y=[y(1) y(2)];
編寫非線性方程組導數的M文件dfx1.m
  function y=dfx1(x)
     y(1)=2*x(1)-10;
     y(2)=2*x(2);
     y(3)=x(2)*x(2)+1;
     y(4)=2*x(1)*x(2)-10;
     y=[y(1) y(2);y(3) y(4)];
     NewtonIterate([0 0])