c
c     Sample_Code_15
c     Solves a nonlinear unsteady state elliptic PDE using Chebyshev Collocation Method
c 
      IMPLICIT DOUBLE PRECISION (a-h,o-z) 
      parameter (NX=16,NY=16,NX1=NX+1,NY1=NY+1)
      DOUBLE PRECISION A(NX1*NY1,NX1*NY1),B(NX1*NY1)
      DOUBLE PRECISION U_G(NX1*NY1),UOLD(NX1*NY1),U(NX1*NY1)
      INTEGER IPVT(NX1*NY1)
      COMMON /CHEB/ X1,X2,Y1,Y2,PI,Alpha,Beta
      COMMON /TRANS/ X(NX1),Y(NY1)
      COMMON /COLLEC/ DX(NX1,NX1),DDX(NX1,NX1),DY(NY1,NY1),DDY(NY1,NY1)
      COMMON /PARAM/ h,BI
      OPEN(11,FILE='sample_code_15_output.txt')
c
      PI    = 4.d0*ATAN(1.d0)
      x1    = 0.d0
      x2    = 1.d0
      y1    = 0.d0
      y2    = 1.d0
      itmax = 10			! Maximum Number of Iterations
      tol   = 1.d-10		! Maximum Tolerance
      iflag = 1			! Flag for the Jacobian Matrix (0 -> analytical Jacobian, 1 -> Numerical Jacobian)
      h     = 1.d-1
      TMAX  = 10.d0
      Time  = 0.d0
c
      call Chebyshev
c
c    Provide Initial Solution
c
      do i = 1, NX1*NY1
          UOLD(i) = 1.d0
      enddo
c
c     Print the Initial solution
c
      write(11,'(A,f10.4)') 'Time =',Time
      write(11,*)
      write(11,'(20x,17f10.4)') (X(i), i = 1, NX1)
      write(11,*)
      do j = 1, NY1
      write(11,'(f10.4,10x,17f10.4)') 
     ! Y(j),(UOLD((i-1)*NY1+j), i = 1, NX1)
      enddo
c
c     Start Integration Loop
c
      count = 0.d0
      istep = 0
  10  Time  = Time + h
      istep = istep + 1
      write(*,*)
      write(*,'(3x,A,f10.6)') 'Time =',Time
      write(*,*)
c
c     Set Termination Criteria for the Integration Loop
      if((Time-TMAX) .gt. 1.d-14) stop
c
c     Provide Inital Guess
c
      do i = 1, NX1*NY1
         U_G(i) = UOLD(i)
      enddo
c
c     Solve the nonlinear PDE iteratively
c
      call Newton(itmax,tol,iflag,U_G,UOLD,U)
c
c     print two dimensional solution
c
      write(11,*)
      write(11,'(3x,A,f10.6)') 'Time =',Time
      write(11,*)
      write(11,'(20x,17f10.4)') (X(i), i = 1, NX1)
      write(11,*)
      do j = 1, NY1
      write(11,'(f10.4,10x,17f10.4)') 
     ! Y(j),(U((i-1)*NY1+j), i = 1, NX1)
      enddo
c
      do i = 1, NX1*NY1
        UOLD(i) = U(i)
      enddo
c
      go to 10
c
1000  stop 
      end
c
c
c
c
c
c
c
c     Subroutine Newton to to solve a system of algebraic equations iteratively
      subroutine Newton(itmax,tol,iflag,UNKNON_G,UOLD,UNKNON)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      parameter (NX=16,NY=16,NX1=NX+1,NY1=NY+1)
      DOUBLE PRECISION A(NX1*NY1,NX1*NY1),B(NX1*NY1)
      DOUBLE PRECISION UNKNON_G(NX1*NY1),UNKNON(NX1*NY1),UOLD(NX1*NY1)
      INTEGER IPVT(NX1*NY1)
      COMMON /CHEB/ X1,X2,Y1,Y2,PI,Alpha,Beta
      COMMON /TRANS/ X(NX1),Y(NY1)
      COMMON /COLLEC/ DX(NX1,NX1),DDX(NX1,NX1),DY(NY1,NY1),DDY(NY1,NY1)

c
c     initialize matrix A and vector B
      do i = 1, NX1*NY1
      UNKNON(i) = UNKNON_G(i)
      B(i) = 0.d0
	do j = 1, NX1*NY1
         A(i,j)=0.d0
	enddo
	enddo
c
c     Start Iterative Solution
      iter = 0
  10  iter = iter + 1
      if (iter .gt. itmax) go to 12
c
c     Input system of nonlinear equations
c     on return B contains the residuals
      call residuals(UNKNON,UOLD,B)
c
c     Check the norm of the residuals
      rnorm = 0.d0
      do i = 1, NX1*NY1
      rnorm = rnorm + B(i)**2
      enddo
      rnorm = dsqrt(rnorm)
      write(*,'(3x,a,i3,5x,a,e17.10)') 
     !                    'iteration #',iter, 'Tolerance =', rnorm
c
      if (rnorm .lt. tol) return
c
c     Evaluate the Jacobian Matrix
      call jacobian(IFLAG,UNKNON,UOLD,A)
c
c     LU-decompose matrix A 
c     Note that matrix A will contains on return the L and U matrixes (original is destroyed)
      call ludcmp(A,NX1*NY1,NX1*NY1,ipvt,dummy)
c
c     Solve the equations: A .(Xiter+1-Xiter) =-B, by backsubstitution
c     Note that vector B will contains the unknowns X (original is destroyed)
c     for Newton Raphson method B = Xiter+1 - Xiter
      do i = 1, NX1*NY1
        B(i) = -B(i)
      enddo
c
      call lubksb(A,NX1*NY1,NX1*NY1,ipvt,B)
c
      do i= 1, NX1*NY1 
        UNKNON(i) = UNKNON(i)+B(i)
      enddo
c
      go to 10
c
  12  write(11,*) 'No Convergence !'
      write(*,*)  'No Convergence !'
      stop
      return
      end
c
c
c
c
c
c
c
c
c
c
c
c
c     Subroutine Residuals to Input System of Non-Linear Equations
      SUBROUTINE RESIDUALS(U,UOLD,R)
      IMPLICIT DOUBLE PRECISION (a-h,o-z) 
      parameter (NX=16,NY=16,NX1=NX+1,NY1=NY+1)
      DOUBLE PRECISION R(NX1*NY1),U(NX1*NY1),UOLD(NX1*NY1)
      COMMON /CHEB/ X1,X2,Y1,Y2,PI,Alpha,Beta
      COMMON /TRANS/ X(NX1),Y(NY1)
      COMMON /COLLEC/ DX(NX1,NX1),DDX(NX1,NX1),DY(NY1,NY1),DDY(NY1,NY1)
      COMMON /PARAM/ h,BI
c
      do i = 1, NX1*NY1
        R(i) = 0.d0
      enddo
c
c     Construct Residuals using Chebyshev Collocation
c
	do i = 2, NX1-1
        do j = 2, NY1-1
          R((i-1)*NY1+j) = R((i-1)*NY1+j) - 
     !    (U((i-1)*NY1+j)-UOLD((i-1)*NY1+j))/h
          do n = 1, NX1
            R((i-1)*NY1+j) = R((i-1)*NY1+j) + DDX(i,n)*U((n-1)*NY1+j)
          enddo
          do m = 1, NY1
            R((i-1)*NY1+j) = R((i-1)*NY1+j) + DDY(j,m)*U((i-1)*NY1+m)
          enddo
        enddo
      enddo
c
c     Implement boundary conditions
c
	do j = 1, NY1 
         R((  1-1)*NY1+j) = U((  1-1)*NY1+j) - 1.d0
         R((NX1-1)*NY1+j) = U((NX1-1)*NY1+j) - 1.d0
      enddo
c
	do i = 1, NX1 
         R((i-1)*NY1+  1) = U((i-1)*NY1+  1) - 1.d0
      enddo
c
	do i = 2, NX1-1 
         R((i-1)*NY1+NY1) = U((i-1)*NY1+NY1) - 0.d0
      enddo
c
      return   
      end 
c
c
c
      SUBROUTINE Chebyshev
      IMPLICIT DOUBLE PRECISION (a-h,o-z) 
      parameter (NX=16,NY=16,NX1=NX+1,NY1=NY+1)
      DOUBLE PRECISION XCHEB(NX1),CX(NX1),CBARX(NX1)
      DOUBLE PRECISION YCHEB(NY1),CY(NY1),CBARY(NY1)
      COMMON /CHEB/ X1,X2,Y1,Y2,PI,Alpha,Beta
      COMMON /TRANS/ X(NX1),Y(NY1)
      COMMON /COLLEC/ DX(NX1,NX1),DDX(NX1,NX1),DY(NY1,NY1),DDY(NY1,NY1)
c
c     Vector XCHEB (1,-1) contains the Gauss-Lobatto grid. 
      do i = 1, NX1 
        XCHEB(i) = dcos((i-1)*pi/NX)
        CBARX(i)  = 1.d0  
        CX(i)     = 1.d0  
      enddo
      do i = 1, NY1 
        YCHEB(i)  = dcos((i-1)*pi/NY)
        CBARY(i)  = 1.d0  
        CY(i)     = 1.d0  
      enddo
c
      CX(1)      = 2.d0  
      CBARX(1)   = 2.d0  
      CBARX(NX1) = 2.d0  
      CY(1)      = 2.d0  
      CBARY(1)   = 2.d0  
      CBARY(NY1) = 2.d0  
c
c     Vector X maps XCHEB to the interval (x1,x2)	==> X(i) = Alpha + Beta*XCHEB(i)
      Alpha1 = (x1+x2)/2.d0
      Beta1  = (x1-x2)/2.d0
      Alpha2 = (y1+y2)/2.d0
      Beta2  = (y1-y2)/2.d0
c
      do i = 1, NX1 
        X(i) = Alpha1 + Beta1*XCHEB(i)
      enddo
c
      do i = 1, NY1 
        Y(i) = Alpha2 + Beta2*YCHEB(i)
      enddo
c
      do i = 1, NX1
         ii = i-1
         do j = 1, NX1
            jj = j-1
            if (ii .ne. jj) then
              DX(i,j) = CBARX(i)/CBARX(j)*(-1)**(ii+jj)
     !                 /(XCHEB(i)-XCHEB(j))
            else
              if (ii .ne. 0 .and. ii .ne. NX) then
              DX(i,j) = -XCHEB(j)/2.d0/(1.d0-XCHEB(j)**2)
              endif
            endif
         enddo
      enddo
      DX( 1, 1) = (2.d0*dble(NX)**2+1.d0)/6.d0
      DX(NX1,NX1) =-(2.d0*dble(NX)**2+1.d0)/6.d0
c
      do i=1,NX1
        do j=1,NX1 
          DDX(i,j) = 0.d0
          do k=1,NX1
            DDX(i,j) = DDX(i,j) + DX(i,k)*DX(k,j)
          enddo
        enddo
      enddo
c
      do i = 1, NX1
         do j = 1, NX1
             DX(i,j) =  DX(i,j)/BETA1
            DDX(i,j) = DDX(i,j)/BETA1/BETA1
         enddo
      enddo
c
      do i = 1, NY1
         ii = i-1
         do j = 1, NY1
            jj = j-1
            if (ii .ne. jj) then
              DY(i,j) = CBARY(i)/CBARY(j)*(-1)**(ii+jj)
     !                 /(YCHEB(i)-YCHEB(j))
            else
              if (ii .ne. 0 .and. ii .ne. NY) then
              DY(i,j) = -YCHEB(j)/2.d0/(1.d0-YCHEB(j)**2)
              endif
            endif
         enddo
      enddo
      DY( 1, 1) = (2.d0*dble(NY)**2+1.d0)/6.d0
      DY(NY1,NY1) =-(2.d0*dble(NY)**2+1.d0)/6.d0
c
      do i=1,NY1
        do j=1,NY1 
          DDY(i,j) = 0.d0
          do k=1,NY1
            DDY(i,j) = DDY(i,j) + DY(i,k)*DY(k,j)
          enddo
        enddo
      enddo
c
      do i = 1, NY1
         do j = 1, NY1
             DY(i,j) =  DY(i,j)/BETA2
            DDY(i,j) = DDY(i,j)/BETA2/BETA2
         enddo
      enddo
c
      return   
      end 
c
c
c
c
c
c
c
c
c     Subroutine jacobian to Evaluate the Jacobian Matrix
      subroutine jacobian(IFLAG,UNKNON,UOLD,A)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      parameter (NX=16,NY=16,NX1=NX+1,NY1=NY+1)
      DOUBLE PRECISION A(NX1*NY1,NX1*NY1),R(NX1*NY1),RJ(NX1*NY1)
      DOUBLE PRECISION UNKNON(NX1*NY1),UNKNONJ(NX1*NY1),UOLD(NX1*NY1)
      COMMON /CHEB/ X1,X2,Y1,Y2,PI,Alpha,Beta
      COMMON /TRANS/ X(NX1),Y(NY1)
      COMMON /COLLEC/ DX(NX1,NX1),DDX(NX1,NX1),DY(NY1,NY1),DDY(NY1,NY1)
      parameter (epsl=1.d-5)
c
      do i = 1, NX1*NY1
	do j = 1, NX1*NY1
        A(i,j)=0.d0
	enddo
	enddo
c
      if (IFLAG .gt. 0) go to 10
c
c     Evaluate the jacobian matrix analatically
      return
c
c     Evaluate the jacobian matrix numerically
 10   do i = 1, NX1*NY1
        UNKNONJ(i) = UNKNON(i)
      enddo
c
      call residuals(UNKNON,UOLD,R)
c
      do i = 1, NX1*NY1
        diff  = dmax1(epsl,dabs(epsl*UNKNON(i)))
        UNKNONJ(i) = UNKNON(i)+diff
        call residuals(UNKNONJ,UOLD,RJ)
        do j = 1, NX1*NY1
          A(j,i) = (RJ(j)-R(j))/diff
        enddo
        UNKNONJ(i)=UNKNON(i)
      enddo
      return
      end
c
c
c
c
c
c
c
c
c
c
c
c
c
c
      SUBROUTINE ludcmp(a,n,np,indx,d)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER n,np,indx(n),NMAX
      DOUBLE PRECISION d,a(np,np),TINY
      PARAMETER (NMAX=500,TINY=1.0e-20)
      INTEGER i,imax,j,k
      DOUBLE PRECISION aamax,dum,sum,vv(NMAX)
      d=1.
      do 12 i=1,n
        aamax=0.d0
        do 11 j=1,n
          if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j))
11      continue
        if (aamax.eq.0.d0) 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.d0
        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.d0)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
c
      SUBROUTINE lubksb(a,n,np,indx,b)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER n,np,indx(n)
      DOUBLE PRECISION a(np,np),b(n)
      INTEGER i,ii,j,ll
      DOUBLE PRECISION 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.d0) 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