title: Receiving NaN value in Fortran app
body: <p>I'm developing an Fortran application for numerically solving Boundary Value Problem for second order ODE of the type: -y''+q(x)*y=r(x). In this application I use Gauss-ellimination algorithm to solve the linear system of equations and write the solution in file. But for the solution vector I receive NaN. Why is that happen? Here is some code.</p>
  
  <pre><code>      subroutine gaussian_solve(s, c, error)
  
      double precision, dimension(:,:), intent(in out) :: s
      double precision, dimension(:),   intent(in out) :: c
           integer        :: error
  
      if(error == 0) then
          call back_substitution(s, c)
      end if
  
       end subroutine gaussian_solve
  !=========================================================================================
  
  !================= Subroutine gaussian_ellimination ===============================
        subroutine gaussion_ellimination(s, c, error)
  
      double precision, dimension(:,:), intent(in out) :: s
      double precision, dimension(:),   intent(in out) :: c
      integer,            intent(out)  :: error
  
      real, dimension(size(s, 1)) :: temp_array
      integer, dimension(1)       :: ksave
      integer                     :: i, j, k, n
      real                        :: temp, m
  
      n = size(s, 1)
  
      if(n == 0) then
          error = -1 
          return
      end if
  
      if(n /= size(s, 2)) then
          error = -2 
          return
      end if
  
      if(n /= size(s, 2)) then 
          error = -3 
          return
      end if
  
      error = 0
      do i = 1, n-1
          ksave = maxloc(abs(s(i:n, i)))
          k = ksave(1) + i - 1
          if(s(k, i) == 0) then
              error = -4
              return
          end if
  
          if(k /= i) then
              temp_array = s(i, :)
              s(i, :) = s(k, :)
              s(k, :) = temp_array
              temp = c(i)
              c(i) = c(k)
              c(k) = temp
          end if
  
          do j = i + 1, n
              m = s(j, i)/s(i, i)
              s(j, :) = s(j, :) - m*s(i, :)
              c(j) = c(j) - m*c(i)
          end do
      end do
  
       end subroutine gaussion_ellimination
       !==========================================================================================
  
     !================= Subroutine back_substitution ========================================
       subroutine back_substitution(s, c)
  
      double precision, dimension(:,:), intent(in) :: s
      double precision, dimension(:),   intent(in out) :: c
  
      real    :: w
      integer :: i, j, n
  
      n = size(c)
  
      do i = n, 1, -1
          w = c(i)
          do j = i + 1, n
              w = w - s(i, j)*c(j)
          end do
          c(i) = w/s(i, i)
      end do
  
        end subroutine back_substitution
  </code></pre>
  
  <p>Where s(i, j) is the matrix of coefficients of the system and c(i) is the solution vector.</p>