c
c     Sample_Code_3
c     Solves a System of IVP's
c 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (N = 2)
      DOUBLE PRECISION Y(N),YOLD(N)
      COMMON x1,x2,h,phi1,phi2
      OPEN(11,FILE='sample_code_3_output.txt')
c
      tol    = 1.d-10		! Maximum Tolerance
      itmax  = 20			! Maximum Number of Iterations
      IFLAG  = 1			! Flag for the Jacobian Matrix (0 -> analytical Jacobian, 1 -> Numerical Jacobian)
      x1     = 0.d0		! Starting Time		
      x2     = 10.d0		! Ending Time
      h      = 1.d-2		! Time Step		
      phi1   = 1.d0		! Parameter
      phi2   = 5.d-1		! Parameter
      IFLAG2 = 3			! 1-> Adams Bashforth, 2-> Predictor Corrector, 3-> Adams Moulten, 4-> RKG4
c
c     provide initial conditions
      YOLD(1) = 1.d0
      YOLD(2) = 0.d0
c
      IF (IFLAG2 .eq. 1) GO TO 10 
      IF (IFLAG2 .eq. 2) GO TO 20 
      IF (IFLAG2 .eq. 3) GO TO 30 
      IF (IFLAG2 .eq. 4) GO TO 40 
c
c     Integrate IVP's using 2nd-order Adams-Bashforth Method
  10  call Adams_Bashforth(N,YOLD,Y)
      stop
c
c     Integrate IVP's using 2nd-order Predictor Corrector Method
  20  call Predictor_Corrector(N,YOLD,Y)
      stop
c
c     Integrate IVP's using 2nd-order Adams-Moulten Method
  30  call Adams_Moulten(N,ITMAX,TOL,IFLAG,YOLD,Y)
      stop
c
c     Integrate IVP's using 4th order Runge-Kutta-Gill Method
  40  call RKG4(N,YOLD,Y)
      stop
      end
c
c
c
c     Subroutine Func to Input RHS functions of IVP
c
c
c
c
      subroutine Func(N,Y,FUN)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION Y(N),FUN(N)
      COMMON x1,x2,h,phi1,phi2
c
      do i = 1, N
        FUN(i) = 0.d0
	enddo
c
      FUN(1) =-phi1*Y(1)
      FUN(2) = phi1*Y(1)-phi2*Y(2)*Y(2)
c
      return
      end
c
c
c
c     Subroutine Adams_Bashforth to integrate IVP's
c
c
c
      subroutine Adams_Bashforth(N,YOLD,Y)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION Y(N),YOLD(N),FUNOLD1(N),FUNOLD2(N)
      COMMON x1,x2,h,phi1,phi2
c
      istep = 0
      Time  = x1
c
c     Print the Initial solution to output file
      write(11,*)
      write(11,'(12x,a,13x,a,13x,a)')  'X','Y(1)','Y(2)'
      write(11,'(3f16.6)') Time,YOLD(1),YOLD(2)
c
c     Start Integration Loop
  10  Time  = Time + h
      istep = istep + 1
c
c     Set Termination Criteria for the Integration Loop
      if((Time-x2) .gt. 1.d-14) return
c     Evaluate the RHS Function for previous time
      call Func(N,YOLD,FUNOLD1)
c
c     For the first time step apply 1st-order explicit Euler
      if (istep .eq. 1) then
        do i = 1, N
          Y(i) = YOLD(i)+h*FUNOLD1(i)
        enddo
      else
c
c     Apply 2nd-order AB Formula for Remaining Time Steps
        do i = 1, N
          Y(i) = YOLD(i)+h/2.d0*(3.d0*FUNOLD1(i)-FUNOLD2(i))
        enddo
      endif
c
c     Print the solution for the current time step to output file
      write(11,'(3f16.6)') Time,Y(1),Y(2)
c
c     Store Previous Time Values
      do i = 1, N
        YOLD(i)    = Y(i)
        FUNOLD2(i) = FUNOLD1(i)
      enddo
c
      go to 10
c
      end
c
c
c
c     Subroutine Predictor-Corrector to integrate IVP's
c
c
c
      subroutine Predictor_Corrector(N,YOLD,Y)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION Y(N),YOLD(N),DUMMY(N),FUNOLD1(N),FUNOLD2(N)
      COMMON x1,x2,h,phi1,phi2
c
      istep = 0
      Time  = x1
c
c     Print the Initial solution to output file
      write(11,*)
      write(11,'(12x,a,13x,a,13x,a)')  'X','Y(1)','Y(2)'
      write(11,'(3f16.6)') Time,YOLD(1),YOLD(2)
c
c     Start Integration Loop
  10  Time  = Time + h
      istep = istep + 1
c
c     Set Termination Criteria for the Integration Loop
      if((Time-x2) .gt. 1.d-14) return
c
c     Evaluate the RHS Function at Previous Time
      call Func(N,YOLD,FUNOLD1)
c
c     For the first time step apply 1st-order explicit Euler for predictor
c     and 1st-order implicit Euler for corrector
      if (istep .eq. 1) then
c       predicor
        do i = 1, N
          Y(i) = YOLD(i)+h*FUNOLD1(i)
        enddo
c       evaluate RHS function at predicted variable
        call Func(N,Y,DUMMY)
c       corrector
        do i = 1, N
          Y(i) = YOLD(i)+h*DUMMY(i)
        enddo
      else
c
c     Apply 2nd-order predictor-corrector method for Remaining Time Steps
c       predicor
        do i = 1, N
          Y(i) = YOLD(i)+h/2.d0*(3.d0*FUNOLD1(i)-FUNOLD2(i))
        enddo
c       evaluate RHS function at predicted variable
        call Func(N,Y,DUMMY)
c       corrector
        do i = 1, N
          Y(i) = YOLD(i)+h/2.d0*(FUNOLD1(i)+DUMMY(i))
        enddo
      endif
c
c     Print the solution for the current time step to output file
      write(11,'(3f16.6)') Time,Y(1),Y(2)
c
c     Store Previous Time Values
      do i = 1, N
        YOLD(i)    = Y(i)
        FUNOLD2(i) = FUNOLD1(i)
      enddo
      go to 10
c
      end
c
c
c
c     Subroutine Adams_Moulten to integrate IVP's
c
c
c
      subroutine Adams_Moulten(N,ITMAX,TOL,IFLAG,YOLD,Y)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION Y(N),Y_G(N),FUN(N),FUNOLD(N),YOLD(N)
      COMMON x1,x2,h,phi1,phi2
c
      istep = 0
      Time  = x1
c
c     Print the Initial solution to output file
      write(11,*)
      write(11,'(12x,a,13x,a,13x,a)')  'X','Y(1)','Y(2)'
      write(11,'(3f16.6)') Time,YOLD(1),YOLD(2)
c
c     Start Integration Loop
  10  Time  = Time + h
      istep = istep + 1
      write(*,'(4x,a,f10.4)')	'Time =',Time
c
c     Set Termination Criteria for the Integration Loop
      if((Time-x2) .gt. 1.d-14) return
c
c     Evaluate the RHS Function for previous time
      call Func(N,YOLD,FUNOLD)
c     
c     Initial guess for current integration step is the solution form previous step
      do i = 1, N
        Y_G(i) = YOLD(i)
      enddo
      call Newton(N,ITMAX,TOL,IFLAG,ISTEP,YOLD,FUNOLD,Y_G,Y)
c
c     Evaluate the RHS Function for current time
      call Func(N,Y,FUN)
c
c     Store previous integration values
      do i = 1, N
        YOLD(i)   = Y(i)
        FUNOLD(i) = FUN(i)
      enddo
c
c     Print the solution for the current time step to output file
      write(11,'(3f16.6)') Time,Y(1),Y(2)
c
      go to 10
c
      return
      end
c
c
c
c     Subroutine RKG4 to integrate IVP's using Rung Kutta Gill 4th order method
c
c
c
      subroutine RKG4(N,YOLD,Y)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION Y(N),FUN(N),FUNOLD(N),YOLD(N)
      DOUBLE PRECISION rk1(N),rk2(N),rk3(N),rk4(N),Y_DUMMY(N),F_DUMMY(N)
      COMMON x1,x2,h,phi1,phi2
c
      a     = (dsqrt(2.d0)-1.d0)/2.d0
      b     = (2.d0-dsqrt(2.d0))/2.d0
      c     = -dsqrt(2.d0)/2.d0
      d     = 1.d0+dsqrt(2.d0)/2.d0
c
      istep = 0
      Time  = x1
c
c     Print the Initial solution to ivp.txt
      write(11,*)
      write(11,'(12x,a,13x,a,13x,a)')  'X','Y(1)','Y(2)'
      write(11,'(3f16.6)') Time,YOLD(1),YOLD(2)
c
c     Evaluate the RHS Function for the Initial Solution
      call Func(N,YOLD,FUNOLD)
c
c     Start Integration Loop
  10  Time  = Time + h
      istep = istep + 1
c
c     Set Termination Criteria for the Integration Loop
      if((Time-x2) .gt. 1.d-14) return
c
      do i = 1, N
        rk1(i)     = h*FUNOLD(i)
        Y_DUMMY(i) = YOLD(i)+0.5d0*rk1(i)
      enddo
      call Func(N,Y_DUMMY,F_DUMMY)
c
      do i = 1, N
        rk2(i)     = h*F_DUMMY(i)
        Y_DUMMY(i) = YOLD(i)+a*rk1(i)+b*rk2(i)
      enddo
      call Func(N,Y_DUMMY,F_DUMMY)
c
      do i = 1, N
        rk3(i)     = h*F_DUMMY(i)
        Y_DUMMY(i) = YOLD(i)+c*rk2(i)+d*rk3(i)
      enddo
      call Func(N,Y_DUMMY,F_DUMMY)
c
      do i = 1, N
        rk4(i)     = h*F_DUMMY(i)
      enddo
c
      do i = 1, N
        Y(i)     = YOLD(i)+1.d0/6.d0*(rk1(i)+rk4(i))+
     !                     1.d0/3.d0*(b*rk2(i)+d*rk3(i))
      enddo
c
      call Func(N,Y,FUN)
c
c     Store previous integration values
      do i = 1, N
        YOLD(i)   = Y(i)
        FUNOLD(i) = FUN(i)
      enddo
c
c     Print the solution for the current time step to ivp.txt
      write(11,'(3f16.6)') Time,Y(1),Y(2)
c
      go to 10
c
      return
      end
c
c
c
c
c     Subroutine Newton to to solve a system of algebraic equations iteratively
c
c
c
      subroutine Newton
     !(N,ITMAX,TOL,IFLAG,ISTEP,YOLD,FUNOLD,UNKNON_G,UNKNON)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION A(N,N),UNKNON_G(N),UNKNON(N),B(N)
      DOUBLE PRECISION YOLD(N),FUNOLD(N)
      integer ipvt(N) 
      COMMON x1,x2,h,phi1,phi2
c
c     initialize matrix A and vector B
      do i = 1, N
        UNKNON(i) = UNKNON_G(i)
	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(N,YOLD,FUNOLD,UNKNON,B)
c
c     Check the norm of the residuals
      rnorm = 0.d0
      do i = 1, N
      rnorm = rnorm + B(i)**2
      enddo
      rnorm = dsqrt(rnorm)
      write(*,'(3x,a,i3,5x,a,e20.10)') 
     !                    'iteration #',iter, 'Tolerance =', rnorm
      if (rnorm .lt. tol) return
c      if (istep .gt. 10) go to 20
c
c     Evaluate the Jacobian Matrix
      call jacobian(IFLAG,N,YOLD,FUNOLD,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,N,N,ipvt,d)
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
  20  do i = 1, N
        B(i) = -B(i)
      enddo
c
      call lubksb(A,N,N,ipvt,B)
c
      do i= 1, N 
        UNKNON(i) = UNKNON(i)+B(i)
      enddo
      go to 10
  12  write(11,*) 'No Convergence !'
      write(*,*)  'No Convergence !'
      stop
      return
      end
c
c
c
c     Subroutine Residuals to Input System of Non-Linear Equations
c
c
c
      subroutine Residuals(N,YOLD,FUNOLD,UNKNON,R)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION UNKNON(N),YOLD(N),FUN(N),FUNOLD(N),R(N)
      COMMON x1,x2,h,phi1,phi2
c
      do i = 1, N
      R(i) = 0.d0
	enddo
c
c     Evaluate the RHS Functions of IVP's for UNKNON
      call Func(N,UNKNON,FUN)
c
c     Evaluate the Residual Equations
      do i = 1, N
        R(i) = UNKNON(i)-YOLD(i)-h/2.d0*(FUN(i)+FUNOLD(i))
      enddo
c
      return
      end
c
c
c
c     Subroutine jacobian to Evaluate the Jacobian Matrix
c
c
c
      subroutine jacobian(IFLAG,N,YOLD,FUNOLD,UNKNON,A)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION A(N,N),UNKNON(N),UNKNONJ(N),FUN(N),FUNJ(N)
      DOUBLE PRECISION YOLD(N),FUNOLD(N)
      COMMON x1,x2,h,phi1,phi2
      parameter (epsl=1.d-5)
c
      do i = 1, N
	do j = 1, N
        A(i,j)=0.d0
	enddo
	enddo
c
      if (IFLAG .gt. 0) go to 10
c
c     Evaluate the jacobian matrix numerically
 10   do i = 1, N
        UNKNONJ(i)=UNKNON(i)
      enddo
      call residuals(N,YOLD,FUNOLD,UNKNON,FUN)
      do i = 1, N
        diff  = dmax1(epsl,dabs(epsl*UNKNON(i)))
        UNKNONJ(i) = UNKNON(i)+diff
        call residuals(N,YOLD,FUNOLD,UNKNONJ,FUNJ)
        do j = 1, N
          A(j,i) = (FUNJ(j)-FUN(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