PROGRAM diff_eq

  USE LU

  IMPLICIT NONE

  REAL, ALLOCATABLE :: x(:), f(:), g(:), y(:), mat(:,:)
  REAL :: dx, a, b, usdx2, us2dx
  INTEGER :: N, i
  REAL, PARAMETER :: y0=1.0, ypN=-1.0/4.0
  CHARACTER (LEN=20) :: fname

  PRINT *,"Inserire estremi a,b dell'intervallo:"
  READ *, a, b
  PRINT *,'Inserire numero di intervalli:'
  READ *, N
  PRINT *,'Nome del file di output:'
  READ *,fname

  ALLOCATE(x(N), f(N), g(N), y(N), mat(N,3))

  dx=(b-a)/N
  DO i=1,N
     x(i)=a+i*dx
  ENDDO
  f=1.0/x
  g=f*f*f
  y=f

  usdx2=1.0/(dx*dx)
  us2dx=1.0/(2.0*dx)
  mat(1,1)=0.0
  mat(1,2)=-2.0*usdx2
  mat(1,3)=usdx2+f(1)*us2dx
  DO i=2,N-1
     mat(i,1)=usdx2-f(i)*us2dx
     mat(i,2)=-2.0*usdx2
     mat(i,3)=usdx2+f(i)*us2dx
  ENDDO
  mat(N,1)=2.0*usdx2
  mat(N,2)=-2.0*usdx2
  mat(N,3)=0.0

  g(1)=g(1)-(usdx2-f(1)*us2dx)*y0
  g(N)=g(N)-2.0*dx*(usdx2+f(N)*us2dx)*ypN

  CALL lu_fact(mat)
  CALL lu_solve(mat,g)

  OPEN(FILE=fname, UNIT=16, STATUS='UNKNOWN')
  WRITE(16,*) a, y0, y0, 0.0
  DO i=1,N
     WRITE(16,*) x(i),g(i),y(i),abs(g(i)-y(i))
  ENDDO
  CLOSE(UNIT=16)

END PROGRAM diff_eq
