c
c     Sample_Code_12
c     Solves a nonlinear parabolic PDE using Chebyshev Collocation Method
c 
c 
      IMPLICIT DOUBLE PRECISION (a-h,o-z) 
      parameter (N=10,N1=N+1)
      DOUBLE PRECISION Y(N1),Y_G(N1),YOLD(N1)
      COMMON /CHEB/ z1,z2,PI,Alpha,Beta
      COMMON /TRANS/ C(N1),CBAR(N1),CT(N1,N1),CTINV(N1,N1),Z(N1)
      COMMON /COLLEC/ D(N1,N1),DD(N1,N1)
      COMMON /PARAM/ Phi,Bi
      COMMON /PARAM2/ TMAX,H,IMETHD,ISTEP
      OPEN(11,FILE='sample_code_12_output.txt')
      OPEN(12,FILE='sample_code_12_output2.txt')
c
      PI     = 4.d0*DATAN(1.d0)
      z1     = 0.d0
      z2     = 1.d0
      Phi    = 1.d0
      Bi     = 1.d0
c
      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)
      IMETHD = 2			! 1 for Exp. Euler, 2 for Imp. Euler
      h      = 1.d-3
      TMAX   = 10.d0
      IPRINT = (TMAX/h)/100
c
      call Chebyshev
c
c    Provide Initial Solution
c
      do i = 1, N1
          YOLD(i) = 1.d0
      enddo
c
c     Print the Initial solution
      write(11,*)
      write(11,'(3x,A,f10.6)') 'Time =',Time
      do i = 1, N1
         write(11,'(10x,2f16.6)') Z(i),YOLD(i)
      enddo
c
c     Start Integration Loop
      count = 0.d0
      istep = 0
      Time  = 0.d0
  10  Time  = Time + h
      istep = istep + 1
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, N1
         Y_G(i) = YOLD(i)
      enddo
c
c     Solve resulting system of nonlinear algebraic equations
c
      call Newton(itmax,tol,iflag,Y_G,YOLD,Y)
c
c     Print the solution for the current time
      if (istep  .eq. iprint) then
        write(11,'(3x,A,f10.6)') 'Time =',Time
        do i = 1, N1
           write(11,'(10x,2f16.6)') Z(i),Y(i)
        enddo
          write(12,'(2f16.6)') Time,Y(N1-1)
          write( *,'(2f16.6)') Time,Y(N1-1)
        istep = 0
      endif
c
      do i = 1, N1
        YOLD(i) = Y(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,YOLD,UNKNON)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      parameter (N=10,N1=N+1)
      DOUBLE PRECISION A(N1,N1),UNKNON_G(N1),UNKNON(N1),B(N1),YOLD(N1)
      COMMON /COLLEC/ D(N1,N1),DD(N1,N1)
      COMMON /PARAM/ Phi,Bi
      COMMON /PARAM2/ TMAX,H,IMETHD,ISTEP
      integer ipvt(N1) 
c
c     initialize matrix A and vector B
      do i = 1, N1
      UNKNON(i) = UNKNON_G(i)
      B(i) = 0.d0
	do j = 1, N1
         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(YOLD,UNKNON,B)
c
c     Check the norm of the residuals
      rnorm = 0.d0
      do i = 1, N1
      rnorm = rnorm + B(i)**2
      enddo
      rnorm = dsqrt(rnorm)
c      write(*,'(3x,a,i3,5x,a,e17.10)') 
c     !                    'iteration #',iter, 'Tolerance =', rnorm
c
      if (rnorm .lt. tol) return
c
c     Evaluate the Jacobian Matrix
      call jacobian(IFLAG,YOLD,UNKNON,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,N1,N1,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, N1
        B(i) = -B(i)
      enddo
c
      call lubksb(A,N1,N1,ipvt,B)
c
      do i= 1, N1 
        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(YOLD,Y,R)
      IMPLICIT DOUBLE PRECISION (a-h,o-z) 
      parameter (N=10,N1=N+1)
      DOUBLE PRECISION YOLD(N1),Y(N1),R(N1)
      COMMON /CHEB/ z1,z2,PI,Alpha,Beta
      COMMON /TRANS/ C(N1),CBAR(N1),CT(N1,N1),CTINV(N1,N1),Z(N1)
      COMMON /COLLEC/ D(N1,N1),DD(N1,N1)
      COMMON /PARAM/ Phi,Bi
      COMMON /PARAM2/ TMAX,H,IMETHD,ISTEP
c
      do i = 1, N1
        R(i) = 0.d0
      enddo
c
c
c     Construct Residuals for the PDE using Chebyshev Collocation
c
      If (IMETHD .eq. 1) then
c     Explicit Euler Integration
      do i = 2, N1-1
        R(i) = (Y(i)-YOLD(i))/h
        do j = 1, N1
          R(i) = R(i) - DD(i,j)*YOLD(j)
        enddo
      enddo
      go to 100
      endif
c
      If (IMETHD .eq. 2) then
c     Implicit Euler Integration
      do i = 2, N1-1
        R(i) = (Y(i)-YOLD(i))/h-Yold(i)
        do j = 1, N1
          R(i) = R(i) - DD(i,j)*Y(j)
        enddo
      enddo
      go to 100
      endif
c
c     Boundary Conditions
c     In Collecation, the first equation for the first BC and the last equation for the second BC
 100  R(1) = 0.d0 
      do j = 1, N1
        R( 1) = R( 1) + D(1,j)*Y(j)
      enddo
      R(N1) = Y(N1) - 0.d0
c
      return   
      end 
c
c
c
c
c
c
c
c
c     Subroutine jacobian to Evaluate the Jacobian Matrix
      subroutine jacobian(IFLAG,YOLD,UNKNON,A)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      parameter (N=10,N1=N+1)
      double precision UNKNON(N1),UNKNONJ(N1),YOLD(N1)
      double precision A(N1,N1),R(N1),RJ(N1)
      parameter (epsl=1.d-5)
      COMMON /PARAM/ Phi,Bi
      COMMON /PARAM2/ TMAX,H,IMETHD,ISTEP
c
      do i = 1, N1
	do j = 1, N1
        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, N1
        UNKNONJ(i) = UNKNON(i)
      enddo
c
      call residuals(YOLD,UNKNON,R)
c
      do i = 1, N1
        diff  = dmax1(epsl,dabs(epsl*UNKNON(i)))
        UNKNONJ(i) = UNKNON(i)+diff
        call residuals(YOLD,UNKNONJ,RJ)
        do j = 1, N1
          A(j,i) = (RJ(j)-R(j))/diff
        enddo
        UNKNONJ(i)=UNKNON(i)
      enddo
      return
      end
c
c
c
      SUBROUTINE Chebyshev
      IMPLICIT DOUBLE PRECISION (a-h,o-z) 
      parameter (N=10,N1=N+1)
      DOUBLE PRECISION XCHEB(N1)
      COMMON /CHEB/ z1,z2,PI,Alpha,Beta
      COMMON /TRANS/ C(N1),CBAR(N1),CT(N1,N1),CTINV(N1,N1),Z(N1)
      COMMON /COLLEC/ D(N1,N1),DD(N1,N1)
c
c     Vector XCHEB (1,-1) contains the Gauss-Lobatto grid. 
      do i = 1, N1 
        XCHEB(i) = dcos((i-1)*pi/N)
        CBAR(i)  = 1.d0  
        C(i)     = 1.d0  
      enddo
c
      C(1)     = 2.d0  
      CBAR(1)  = 2.d0  
      CBAR(N1) = 2.d0  
c
c     Vector X maps XCHEB to the interval (z1,z2)	==> X(i) = Alpha + Beta*XCHEB(i)
      Alpha = (z1+z2)/2.d0
      Beta  = (z1-z2)/2.d0
c
      do i = 1, N1 
        Z(i) = Alpha + Beta*XCHEB(i)
      enddo
c
      do i = 1, N1
      do j = 1, N1
        CT(i,j)    = 2.d0/N/CBAR(i)/CBAR(j)*dcos(pi*(i-1)*(j-1)/N)  
        CTINV(i,j) = dcos(pi*(i-1)*(j-1)/N)   
      enddo
      enddo
c
      do i = 1, N1
         ii = i-1
         do j = 1, N1
            jj = j-1
            if (ii .ne. jj) then
              D(i,j) = CBAR(i)/CBAR(j)*(-1)**(ii+jj)/(XCHEB(i)-XCHEB(j))
            else
              if (ii .ne. 0 .and. ii .ne. N) then
              D(i,j) = -XCHEB(j)/2.d0/(1.d0-XCHEB(j)**2)
              endif
            endif
         enddo
      enddo
      D( 1, 1) = (2.d0*dble(N)**2+1.d0)/6.d0
      D(N1,N1) =-(2.d0*dble(N)**2+1.d0)/6.d0
c
      do i=1,N1
        do j=1,N1 
          DD(i,j) = 0.d0
          do k=1,N1
            DD(i,j) = DD(i,j) + D(i,k)*D(k,j)
          enddo
        enddo
      enddo
c
      do i = 1, N1
         do j = 1, N1
             D(i,j) =  D(i,j)/BETA
            DD(i,j) = DD(i,j)/BETA/BETA
         enddo
      enddo
c
      return   
      end 
c
c
c
      SUBROUTINE Derivative(YHAT,D1YHAT,D2YHAT)
      IMPLICIT DOUBLE PRECISION (a-h,o-z) 
      parameter (N=10,N1=N+1)
      DOUBLE PRECISION YHAT(N1),D1YHAT(N1),D2YHAT(N1)
      COMMON /CHEB/ z1,z2,PI,Alpha,Beta
      COMMON /TRANS/ C(N1),CBAR(N1),CT(N1,N1),CTINV(N1,N1),Z(N1)
c
      do i = 1, N1 
        D1YHAT(i) = 0.d0
        D2YHAT(i) = 0.d0
      enddo
c
      do i = 1, N1
        ri = dble(i) - 1.d0
c       1st-oder derivatives
        do ip = i + 1, N1, 2
          rp = dble(ip) - 1.d0
          D1YHAT(i) = D1YHAT(i) + (2.d0/C(i)*rp)/Beta*YHAT(ip)
        enddo
c       2nd-oder derivatives
        do ip = i + 2, N1, 2
          rp = dble(ip) - 1.d0
          D2YHAT(i) = D2YHAT(i) + 
     !    1.d0/C(i)*rp*(rp*rp-ri*ri)/Beta/Beta*YHAT(ip)
        enddo
      enddo
c
      return   
      end 
cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
      subroutine RESP(AR,AS) 
cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
      implicit double precision (a-h,o-z)
      parameter (N=10,N1=N+1)
      DOUBLE PRECISION AR(N1),AS(N1)
      COMMON /TRANS/ C(N1),CBAR(N1),CT(N1,N1),CTINV(N1,N1),Z(N1)
c
c     invert chebychev
c
      do i = 1, N1 
        AS(i) = 0.d0 
        do j = 1, N1 
          AS(i) = AS(i) + CT(i,j) * AR(j)
        enddo
      enddo
c
      return  
      end 
cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
      subroutine SPRE(AS,AR) 
cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
      implicit double precision (a-h,o-z)
      parameter (N=10,N1=N+1)
      DOUBLE PRECISION AR(N1),AS(N1)
      COMMON /TRANS/ C(N1),CBAR(N1),CT(N1,N1),CTINV(N1,N1),Z(N1)
c
c     invert chebychev
c
      do i = 1, N1 
        AR(i) = 0.d0 
        do j = 1, N1 
          AR(i) = AR(i) + CTINV(i,j) * AS(j)
        enddo
      enddo
c
      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