Program M_Newton_Num_Sys   ! main pt. FIXS - "fix" pt. sistem neliniar
real(4), allocatable :: x(:),f(:)
external fcns
character(260) nume, fname
character(80) text(-1:1)
character*1 w
data delta /1E-3/
integer width
character(6) str_n	!, str*11
logical print_it

character(132), allocatable :: f_text(:)
character(132) f_source
character(2) i_char

text(-1)='** Matrice singulara'
text(0) ='* XTOL atins'
text(1) ='** LNIT depasit'

print '(a, \)', ' Introdu numarul de ecuatii: '
read *, n

allocate(f_text(n))
! Initialize with "f(i) = ?":
do i=1,n
iwi =width(i)
write(i_char, '(I<iwi>)') i
f_text(i) ='f('//trim(i_char)//') = ?'
enddo

! Open "f_source" and read:
f_source ='fcns.f90'
call F_Expr_Sys(f_source, n, f_text)

do i=1,n
write(*, '(1x, a)') trim(f_text(i))
enddo

print '(a, \)', ' Nume fisier de iesire (fara extensie): '
read(*,'(a)') nume
fname =trim(nume)//'.rez'
open(1, file= fname, status ='unknown')

allocate (x(n), f(n))

write(1, '(a)') ' Newton-numeric: Jacobianul calculat numeric'

do i=1,n
if(i==1) write(1,*)

write(1, '(1x, a)') trim(f_text(i))
enddo

deallocate(f_text)

nrad =1

iwn =width(n)
write(str_n, '(I<iwn>)') n

!print '(a, \)',' Introdu aproximatia initiala "x(1:n)": '
100 print '(a, \)',' Introdu aproximatia initiala x(1:'//trim(str_n)//'): '
read*,x

print '(a, \)',' Introdu xtol, lnit: '
read*, xtol, lnit

print '(a, \)',' Cod de tiparire=? (0/#0): '
read*, kod
print_it = kod /=0

iwnrad =width(nrad); iwlnit =width(lnit);
write(1, '(/, a, I<iwnrad>, a1)') ' Radacina nr.', nrad

write(1,*) 'Aproximatia initiala:', x
write(1, '(a, 1pG10.3, a, I<iwlnit>)') ' xtol:', xtol, '; lnit: ', lnit


call Newton_Sys_Num(fcns, n, x, xtol, delta, lnit, f, kod)

write(1,'(/, a)') trim(text(kod))

iwlnit =width(lnit)
if(kod ==0) then
write(1,*) ' Solutia X: ', x
write(1,*) ' F(X):      ', f
write(1, '(a, I<iwlnit>)') '  Numar de iteratii = ', lnit

else
write(1, '(a, I<iwlnit>)') '  Iteratia numarul: ', lnit
write(1,*) ' Aproximatia curenta: ', x
write(1,*) ' F(X): ', f
endif


print '(/, a, \)',' Alta radacina ? (Y,y/char): '
read '(a1)', w

if(w.eq.'Y'.or.w.eq.'y') then
nrad =nrad+1
goto 100
endif

deallocate(x,f)
close(1)

print '(/, a, a, /)', ' * Fisier de rezultate: ', trim(fname)

end

