
(1) 

(2) 
> 
f:=array([seq(evalf(2*x[i]+x[i]^2+0.99/i^2),i=1..N)]);


(3) 
> 
fsolve({f[1],f[2],f[3],f[4]},{x[1]=1.2,x[2]=1,x[3]=0,x[4]=0});


(4) 
> 
eqsA := [seq(F[i]=f[i],i=1..N)]:

> 
irform2 := StatSeq(seq(Assign(op(i)),i=eqsA)):

> 
prccons:= codegen[intrep2maple](Proc(Parameters(x::Array,F::Array),irform2)):

> 
f3:=Compiler:Compile(prccons):

> 
x0:=Vector([seq(1/2,i=1..N)],datatype=float[8]);


(5) 
> 
f0:=Vector([seq(0,i=1..N)],datatype=float[8]);


(6) 
> 
eqsJ := [seq(seq(Jac[i,j]=evalf(diff(f[i],x[j])),j=1..N),i=1..N)];


(7) 
> 
irformJ := StatSeq(seq(Assign(op(i)),i=eqsJ)):

> 
prcconsJ:= codegen[intrep2maple](Proc(Parameters(x::Array,Jac::Array),irformJ)):

> 
J3:=Compiler:Compile(prcconsJ):

> 
j0:=Matrix(1..N,1..N,datatype=float[8]):

Following linear sovler algorithm is basiclly from dsolve/numeric/SC
.
> 
s3:=proc(n::posint, A::Array( datatype = float[ 8 ] ) , ip::Array( datatype = integer[ 4 ] ),b::Array( datatype = float[ 8 ] ) ) local i::integer, j::integer, k::integer, m::integer, t::float; ip[n] := 1; for k to n1 do m := k; for i from k+1 to n do if abs(A[m,k]) < abs(A[i,k]) then m := i end if end do; ip[k] := m; if m <> k then ip[n] := ip[n] end if; t := A[m,k]; A[m,k] := A[k,k]; A[k,k] := t; if t = 0 then ip[n] := 0; #return ip end if; for i from k+1 to n do A[i,k] := A[i,k]/t end do; for j from k+1 to n do t := A[m,j]; A[m,j] := A[k,j]; A[k,j] := t; if t <> 0 then for i from k+1 to n do A[i,j] := A[i,j]+A[i,k]*t end do end if end do end do; if A[n,n] = 0 then ip[n] := 0 end if; #ip[n] if ip[n] = 0 then error "The matrix A is numerically singular" end if; if 1 < n then for k to n1 do m := ip[k]; t := b[m]; b[m] := b[k]; b[k] := t; for i from k+1 to n do b[i] := b[i]+A[i,k]*t end do end do; for k from n by 1 to 2 do b[k] := b[k]/A[k,k]; t := b[k]; for i to k1 do b[i] := b[i]+A[i,k]*t end do end do end if; b[1] := b[1]/A[1,1]; return; end proc:

> 
s3c:=Compiler:Compile(s3):

> 
x0:=Vector(N,[seq(0.5,i=1..N)],datatype=float[8]);


(8) 

(9) 
> 
j0:=Matrix(1..N,1..N,datatype=float[8]);


(10) 

(11) 
> 
Newton:=proc(x0,f0,j0,N) global f3,s3c; local xnew,err,xold,ip,i::integer; ip:=Vector(N,[seq(1,i=1..N)],datatype=integer[4]); err:=10; xold:=x0; while err>1e8 do f3(xold,f0); J3(xold,j0); s3c(N,j0,ip,f0); xnew:=xold+f0; err:=max(abs(f0))/N; xold:=xnew; end: xnew; end proc;


(12) 
> 
Compiler:Compile(Newton);

Error, (in CodeGeneration:IssueError) cannot analyze noninteger range boundary args[4]


> 
t11:=time():Newton(x0,f0,j0,N);time()t11;


(13) 
