PROGRAM TP4
  IMPLICIT NONE
  INTEGER :: N, no, in, k, Meth, m
  REAL    :: Dt, x0, Ert, C, Tmax, Xk, Xk1, Xk2, Xk3, Xk4
  REAL, DIMENSION(4) :: ErrL2, ErrLI

  Tmax = 50
  C    = 0.20

  ! Compute the order of accurecy of the scheme
  ! using the exact solution of dX/dt = -c X
  ! -----------------------------------------------------------
  DO Meth = 1, 6
     N = 20
     SELECT CASE(Meth)
        CASE(1) ! Explicit Euler 
        WRITE(6,*) " Explicit Euler scheme "

        CASE(2) ! Runge-Kutta RK2
        WRITE(6,*) " Runge-Kutta RK2 scheme "

        CASE(4) ! Runge-Kutta RK4
        WRITE(6,*) " Runge-Kutta RK4 scheme "

     CASE DEFAULT
        ! Runge-Kutta RKp
        WRITE(6,'(A,I2)') " Runge-Kutta RKp scheme p =",Meth

     END SELECT

     DO no = 1,4,1
        N = N*2
        X0 = 1.0
        Dt = Tmax/N
        ErrL2(no) = 0.0
        ErrLI(no) = 0.0
        SELECT CASE(Meth)
           CASE(1) ! Explicit Euler 
           DO in = 1, N

              X0 = X0 - C*dt*X0

              ErrL2(no) = ErrL2(no) + (  X0- EXP( - c*in*dt ))**2
              ErrLI(no) = MAX(ErrLI(no), ABS( X0- EXP( - c*in*dt )) )
           END DO

           CASE(2) ! Runge-Kutta RK2
           DO in = 1, N

              Xk = X0 - 0.5*C*dt*X0
              X0 = X0 - C*dt*Xk

              ErrL2(no) = ErrL2(no) + (  X0- EXP( - c*in*dt ))**2
              ErrLI(no) = MAX(ErrLI(no), ABS( X0- EXP( - c*in*dt )) )
           END DO

           CASE(4) ! Runge-Kutta RK2
           DO in = 1, N

              Xk2 = X0 - 0.5*C*dt*X0
              Xk3 = X0 - 0.5*C*dt*Xk2
              Xk4 = X0 -     C*dt*Xk3

              X0  = X0 - C*dt*( X0 + 2.0*Xk2 + 2.0*Xk3 + Xk4)/6.0

              ErrL2(no) = ErrL2(no) + (  X0- EXP( - c*in*dt ))**2
              ErrLI(no) = MAX(ErrLI(no), ABS( X0- EXP( - c*in*dt )) )
           END DO
        CASE DEFAULT

           DO in = 1, N
              Xk = X0
              DO m = 0, Meth-1
                 Xk = X0 - C*dt*Xk/( Meth - m)
              END DO


              X0  = Xk

              ErrL2(no) = ErrL2(no) + (  X0- EXP( - c*in*dt ))**2
              ErrLI(no) = MAX(ErrLI(no), ABS( X0- EXP( - c*in*dt )) )
           END DO

        END SELECT

        ErrL2(no) = SQRT(ErrL2(no)/N)
        WRITE(6, '( 10("-"), "  N=",I5 )' ) N

     END DO



     WRITE(6,*) " Estimated L2     Order of accuracy = ", LOG( ErrL2(1)/ErrL2(2))/LOG(2.0) 
     WRITE(6,*) " Estimated L2     Order of accuracy = ", LOG( ErrL2(2)/ErrL2(3))/LOG(2.0) 
     WRITE(6,*) " Estimated L2     Order of accuracy = ", LOG( ErrL2(3)/ErrL2(4))/LOG(2.0) 

     WRITE(6, '( 30("-") )' ) 

     WRITE(6,*) " Estimated Linfty Order of accuracy = ", LOG( ErrLI(1)/ErrLI(2))/LOG(2.0) 
     WRITE(6,*) " Estimated Linfty Order of accuracy = ", LOG( ErrLI(2)/ErrLI(3))/LOG(2.0) 
     WRITE(6,*) " Estimated Linfty Order of accuracy = ", LOG( ErrLI(3)/ErrLI(4))/LOG(2.0) 

     WRITE(6, '( 70("-") )' ) 
  END DO


END PROGRAM TP4
