Метод Ньютона для нелинейной системы уравнений

Я новичок в программировании на Fortran90. Я использовал метод NR для системы нелинейных уравнений, найденной в Numerical Recipes, и собрал код, который не генерирует ошибок при компиляции с помощью GFortran. Проблема в том, что он также не генерирует никакого значения для вывода. Может ли это быть из-за того, что мое первоначальное предполагаемое корневое значение далеко от фактического корня, или я сделал ошибку в этом коде?

Любая помощь/совет по этому вопросу будет высоко оценена.

  program main
    implicit real*8(a-h,o-z)
    parameter(n=4)   
    logical check
    real*8 x(n),fvec(n),fjac(n)
    open(20,file="output1.txt",status="unknown")
    do i=1,n
      x(i)= 4.
    enddo
    call mnewt(ntrial,x,n,tolx,tolf)    
        call usrfun(x,n,fvec,fjac)
    do i=1,n
      write(20,*) 'x(',i,')=',x(i),fvec(i)
    enddo
  end

subroutine mnewt(ntrial,x,n,tolx,tolf)

integer n,ntrial,np
real*8 tolf,tolx,x(n)
parameter (np=15)
!uses lubksb, ludcmp, usrfun
! Given an initial guess x for a root in n dimensions, take ntrial Newton-Raphson steps to
! improve the root. Stop if the root converges in either summed absolute variable increments
! tolx or summed absolute function values tolf.
integer i,k,indx(np)
real*8 d,errf,errx,fjac(np,np),fvec(np),p(np)

do 14 k=1,ntrial
  call usrfun(x,n,fvec,fjac)
  errf=0.
  do 11 i=1,n
      errf=errf+abs(fvec(i))
  11 continue
  if(errf.le.tolf)return
  do 12 i=1,n
      p(i)=-fvec(i)
  12 continue

  call ludcmp(fjac,n,np,indx,d)
  call lubksb(fjac,n,np,indx,p)

  errx=0.
  do 13 i=1,n
     errx=errx+abs(p(i))
     x(i)=x(i)+p(i)
   13 continue
   if(errx.le.tolx)return

14 continue
return
end

subroutine usrfun(x,n,fvec,fjac)

    implicit none

    integer n
    real*8 x(n),fvec(n),fjac(n,n), hl, ul, br, bl   

hl=1.00
ul=1.00
br=0.20
bl=0.00

! Initial guesses

        x(1)=0.0
        x(2)=1.5
        x(3)=0.5
        x(4)=0.5

    fvec(1)=(x(2))+(2*sqrt((x(1))))-ul-(2*(sqrt(hl)))
    fvec(2)=((x(3))*(x(4)))-((x(1))*(x(2)))
    fvec(3)=((x(3))*(x(4))*(x(4)))+(0.5*(x(3))*(x(3)))-((x(1))*(x(2))*(x(2)))-(0.5*(x(1))*(x(1)))+(0.5*(br-bl)*x(1)+x(3))
    fvec(4)=(x(4))-sqrt((x(3)))

    fjac(1,1)=((x(1))**(-0.5))
    fjac(1,2)=1
    fjac(1,3)=0
    fjac(1,4)=0
    fjac(2,1)=(-x(2))
    fjac(2,2)=(-x(1))
    fjac(2,3)=x(4)
    fjac(2,4)=x(3)
    fjac(3,1)=((x(2))**2)-(x(1))+(0.5)*(br-bl)
    fjac(3,2)=-2*((x(1))*(x(2)))
    fjac(3,3)=((x(4))*(x(4)))+(x(3))+(0.5)*(br-bl)*(x(3))
    fjac(3,4)=2*((x(3))*(x(4)))
    fjac(4,1)=0
    fjac(4,2)=0
    fjac(4,3)=-0.5*((x(3))**(-0.5))
    fjac(4,4)=1

end subroutine usrfun


     subroutine ludcmp(a,n,np,indx,d)  !fjac=a
      integer n,np,indx(n),nmax
      real*8 d,a(np,np),tiny
      parameter (nmax=2500,tiny=1.0e-20)
      integer i,imax,j,k
      real*8 aamax,dum,sum,vv(nmax)

      d=1.
      do 12 i=1,n
        aamax=0.
        do 11 j=1,n
          if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j))
!   print*,a(21,j)
11      continue
!      print*,i,aamax
!   pause
        if (aamax.eq.0.) pause 'singular matrix in ludcmp'
        vv(i)=1./aamax
12    continue
      do 19 j=1,n
        do 14 i=1,j-1
          sum=a(i,j)
          do 13 k=1,i-1
            sum=sum-a(i,k)*a(k,j)
13        continue
          a(i,j)=sum
14      continue
        aamax=0.
        do 16 i=j,n

          sum=a(i,j)
          do 15 k=1,j-1
            sum=sum-a(i,k)*a(k,j)
15        continue
          a(i,j)=sum
          dum=vv(i)*abs(sum)
          if (dum.ge.aamax) then
            imax=i
            aamax=dum
          endif
16      continue
        if (j.ne.imax)then
          do 17 k=1,n
            dum=a(imax,k)
            a(imax,k)=a(j,k)
            a(j,k)=dum
17        continue
          d=-d
          vv(imax)=vv(j)
        endif
        indx(j)=imax
        if(a(j,j).eq.0.)a(j,j)=tiny
        if(j.ne.n)then
          dum=1./a(j,j)

          do 18 i=j+1,n
            a(i,j)=a(i,j)*dum
18        continue
        endif
19    continue
      return
      end

!lubksb

     subroutine lubksb(a,n,np,indx,b)
      integer n,np,indx(n)
      real*8 a(np,np),b(n)
      integer i,ii,j,ll
      real*8 sum
      ii=0
      do 12 i=1,n
        ll=indx(i)
        sum=b(ll)
        b(ll)=b(i)
        if (ii.ne.0)then
          do 11 j=ii,i-1
            sum=sum-a(i,j)*b(j)
11        continue
        else if (sum.ne.0.) then
          ii=i
        endif
        b(i)=sum
12    continue
      do 14 i=n,1,-1
        sum=b(i)
        do 13 j=i+1,n
          sum=sum-a(i,j)*b(j)
13      continue
        b(i)=sum/a(i,i)
14    continue
      return
      end