c=================================
      program main
      implicit none
      integer n
      parameter (n=5)
      double precision a(n,n), f(n), u(n)
      double precision h
      integer i,j
      integer pivot(n),ok
      
      real myfunction
c================================
      h = 1D0/(n+1)
c just an example how to set up calling a function defined in 
c a different file/module
      
      call myroutine()
      write (6,*) myfunction()

c better safe than sorry: initialize crucial variables
      do i = 1, n
         do j=1, n
            a(i,j) = 0D0
         enddo
         f(i) = 0D0
      enddo

c set up the matrix and right hand side
      do i = 1, n
         a(i,i) = 2D0
         if (i.gt.1) a(i,i-1) = -1D0
         if (i.lt.n) a(i,i+1) = -1D0
         f (i) = h*h;
      enddo
c *** can you think of how to improve the loop above ? 
c *** consider the cost of an "if" statement as equal to a flop

c call the linear solver from lapack family www.netlib.org
      call DGESV(n, 1, a, n, pivot, f, n, ok)

c it is a good idea to see if the solver completed
      write(6,*) ok

c solution is returned in the vector f
      write(6,*) (f(i),i=1,n)

c *** what else is returned ? check in the description of the routine
c *** see the factors that a(.,.) was factored to 

c *** i) can you find a solver better suited to this problem ?
c ***    hint: this matrix is ...  and ... 

c *** PUSH the code to solve the problem as large as possible

      end

c -----------------------------------
      real function myfunction()
c *** can you fix this code ? it is clearly not giving the right value !
      myfunction = 1/3
      end