PROGRAM QUESTION3
INTEGER I,NI,J,NJ,N,NP,NMAX,II,IL,IR
PARAMETER (NI=10,NJ=10)
REAL DX,DY,XTOTAL,YTOTAL,D,TINY, dt, tfin, t
PARAMETER(NP=100,N=100,NMAX=100)
DIMENSION A(NP,NP),INDEX(N),B(N)
REAL XNI,YNJ

        XNI=10.0 
        YNJ=10.0
        
        XTOTAL=1.
        YTOTAL=1.
              
        DX=XTOTAL/(XNI-1)
        DY=YTOTAL/(YNJ-1)
		

!INITIALISE MATRIX AND RHS

DO I=1,N
        DO J=1,N
           A(I,J)=0.
        ENDDO
        B(I)=0.
ENDDO  
!COMPLETED

!FILL IN MATRIX WITH COEFFICIENTS FOR INTERNAL NODES

DO I=2,NI-1
        DO J=2,NJ-1
                IL=NI*(J-2)+I
                II=NI*(J-1)+I
                IR=NI*(J)+I
                A(II,IL)=1./(DY*DY)
                A(II,II-1)=1./(DX*DX)
                A(II,II)=-2.*(1./(DX*DX)+1./(DY*DY))
                A(II,II+1)=1./(DX*DX)
                A(II,IR)= 1./(DY*DY)
        ENDDO
ENDDO
!COMPLETED  

 ! LEFT SIDE
   I=1
DO J=1,NJ-1
        II=NI*(J-1)+I
        A(II,II)=1.
        B(II)=300
ENDDO
! CORNER NODE 
        I=1
        J=NJ
        II=NI*(J-1)+I
        A(II,II)=1.
        B(II)=350
!COMPLETED

!TOP SIDE
J=NJ
DO I=2,NI
        II=NI*(J-1)+I
        A(II,II)=1.
        B(II)=400
ENDDO
!COMPLETED

!BOTTOM SIDE                       
   J=2
DO I=2,NI-1
        II=I
        IR=NI*(J-1)+I
        A(II,II)=1./(DY)
        A(II,IR)=-1./(DY)
ENDDO
!COMPLETED

!CORNER NODE
A(10,10)=1. ! modify here
A(10,19)=-1. ! modify here
!COMPLETED

!RIGHT SIDE                
DO J=2,NJ-1
        II=NI*(J-1)+NI
        IL=NI*(J-1)+(NI-1)
        A(II,IL)=-1./(DX)
        A(II,II)=1./(DX)
ENDDO
!COMPLETED 

                   
  !INITIALISE MATRIX AND RHS
       
        CALL LUDCMP(A,N,NP,INDEX,D)
        
        CALL LUBKSB(A,N,NP,INDEX,B)    
    
        DO I=1,N
                WRITE(*,*) B(I)
        ENDDO
        
         OPEN(1,FILE='results.txt') 
         DO I=1,NI
                 DO J=1,NJ
               II=NI*(J-1)+I
                 WRITE(1,*)B(II)
                 ENDDO
         ENDDO
        
END PROGRAM               
      
!------------------------------------------------

SUBROUTINE LUDCMP( A, N, NP, INDEX, D )

  INTEGER N, NP, INDEX(N), NMAX
  REAL D,A(NP,NP),TINY
  PARAMETER (NMAX=100,TINY=1.0E-20)
  INTEGER I,IMAX,J,K
  REAL  AAMAX,DUM,SUM,VV(NMAX)

   D=1.

   DO I=1,N
     AAMAX=0.
     DO J=1,N

     IF (ABS(A(I,J)).GT.AAMAX)AAMAX=ABS(A(I,J))

     ENDDO

     IF (AAMAX.EQ.0.) PAUSE 'SINGULAR MATRIX ON LUDCOMP'

     VV(I)=1./AAMAX

  ENDDO

  DO J=1,N
    DO I=1,J-1
       SUM=A(I,J)

       DO K=1,I-1
         SUM=SUM-A(I,K)*A(K,J)
       ENDDO

      A(I,J)=SUM
    ENDDO

    AAMAX=0.

    DO I=J,N
      SUM=A(I,J)

      DO K=1,J-1
        SUM=SUM-A(I,K)*A(K,J)
      ENDDO

      A(I,J)=SUM

      DUM=VV(I)*ABS(SUM)

      IF(DUM.GE.AAMAX)THEN
            IMAX=I
            AAMAX=DUM
      ENDIF
    ENDDO

    IF (J.NE.IMAX) THEN
      DO K=1,N
              DUM=A(IMAX,K)
              A(IMAX,K)=A(J,K)
              A(J,K)=DUM
      ENDDO
      D=-D
      VV(IMAX)=VV(J)
    ENDIF

    INDEX(J)=IMAX

    IF(A(J,J).EQ.0.)A(J,J)=TINY

    IF(J.NE.N)THEN
            DUM=1./A(J,J)
            DO I=J+1,N
            A(I,J)=A(I,J)*DUM
            ENDDO
    ENDIF

  ENDDO

  RETURN
END
!------------------------------------------------
!------------------------------------------------
SUBROUTINE LUBKSB( A, N, NP, INDEX, B )

  INTEGER N,NP,INDEX(N)
  REAL  A(NP,NP),B(N)
  INTEGER I,II,J,LL
  REAL  SUM

  II=0

  DO I=1,N
  LL=INDEX(I)
  SUM=B(LL)
  B(LL)=B(I)

      IF(II.NE.0.)THEN
        DO J=II,I-1
          SUM=SUM-A(I,J)*B(J)
        ENDDO
      ELSEIF(SUM.NE.0.) THEN
         II=I
      END IF

  B(I)=SUM
  ENDDO

  DO I = N, 1, -1
  SUM=B(I)

  DO J=I+1,N
          SUM=SUM-A(I,J)*B(J)
  END DO
  B(I)=SUM/A(I,I)
  ENDDO

  RETURN
END
!------------------------------------------------