title: Newton Method for Nonlinear set of Equations
body: <p>I'm a rookie in programming in Fortran90. I used NR method for a system of non-linear equations found in Numerical Recipes and put together a code that does not generate any errors when I compile with GFortran. Problem is, it does not generate any value for output either.
  Could it be because my initial guess root value far off from actual root or have I made an error in this code?</p>
  
  <p>Any help/advice on this matter will be highly appreciated.</p>
  
  <pre><code>  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
  </code></pre>