-- This pad solves Newton's Method for a system of non-linear equations.
--by Anthony Egan, M.S.
--SD School of Mines
--afe8607@silver.sdsmt.edu

-- *** Requires SVD XFun for matrix inversion ***

multVector(A,vec)[i] = sum(A[i,j]*vec[j],j,1,count(vec)) dim[count(A)]
zero2[i,i]=0 dim[n,n]

--- 
n=3       -- number of equations

f={10*x[1]^3-5*x[1]-x[2],
   5*x[2]^3-3*x[2]-x[3],
   2*x[3]^3-x[3]-8}  -- function matrix

x:={1,1,1}:-- guess of vector x

tol=1e-5  -- average difference between final iterations


df:=zero2:
(i:=1,
 (j:=1,
  (x[j]:=1.01*x[j],f2:=f,
   x[j]:=x[j]/1.01,f1:=f,
   df[i,j]:=(f2[i]-f1[i])/0.01/x[j],
   j:=j+1) while j<=n,
 i:=i+1) while i<=n,
 invert(df),,
 x:=x-multVector(_Inv,f)) while sum(f[at],at,1,n)^2 >= tol*n:

f:{0.000,0.000,0.000}--should be close to zero
x:{0.789,0.973,1.692}--final solution