!********************************************************************
!*    Solve Y' = F(X,Y) with Initial Condition Y(X0)=Y0 using the   *
!*    Euler-Romberg Method                                          *
!* ---------------------------------------------------------------- *
!* REFERENCE: "Méthode de calcul numérique- Tome 2 - Programmes en  *
!*             Basic et en Pascal By Claude Nowakowski, Edition du  *
!*             P.S.I., 1984" [BIBLI 04].                            * 
!* ---------------------------------------------------------------- *
!* SAMPLE RUN:                                                      *
!* (Solve Y' = X*X*Y with Y(0)=1 for X0=0 up to X1=1.0, exact       *
!*  solution is Y = Exp(X^3/3) ).                                   *
!*                                                                  *
!* Starting X...: 0                                                 *
!* Ending X.....: 1                                                 *
!* Value Y at X0: 1                                                 *
!* Initial step : 0.1                                               *
!* Desired error: 1e-6                                              *
!*                                                                  *
!*    X         Y            Y       Error       Number of          *
!*          estimated      exact                subdivisions        *
!* ----------------------------------------------------------       *
!*   0.1     1.000333    1.000333  0.00000000       4               *
!*   0.2     1.002670    1.002670  0.00000001       4               *
!*   0.3     1.009041    1.009041  0.00000006       4               *
!*   0.4     1.021562    1.021562  0.00000014       4               *
!*   0.5     1.042547    1.042547  0.00000027       4               *
!*   0.6     1.074654    1.074655  0.00000086       4               *
!*   0.7     1.121125    1.121126  0.00000107       4               *
!*   0.8     1.186094    1.186095  0.00000126       4               *
!*   0.9     1.275067    1.275069  0.00000133       4               *
!*   1.0     1.395611    1.395612  0.00000114       4               *
!*   1.1     1.558410    1.558411  0.00000047       4               *
!* ----------------------------------------------------------       *
!*                                                                  * 
!*                                F90 Release By J-P Moreau, Paris. *
!*                                       (www.jpmoreau.fr)          *
!********************************************************************
Program EulerRomberg

implicit real*8 a-h,o-z

parameter(NMAX=100)

  DIMENSION X(NMAX), Y(NMAX), NL(NMAX)

  !Initial data
  print *,' '
  write(*,10,advance='no'); read *, X0
  write(*,20,advance='no'); read *, X1
  write(*,30,advance='no'); read *, Y0
  write(*,40,advance='no'); read *, H
  write(*,50,advance='no'); read *, ER

  !number of calculations
  NC=Int((X1-X0)/H)

  !call integration routine
  call Euler_Romberg(NC,H,ER,X0,Y0,X,Y,NL)

  !write header
  PRINT *,' '
  PRINT *,'   X         Y            Y       Error      Number of   '
  PRINT *,'         estimated      exact               subdivisions '
  PRINT *,'---------------------------------------------------------'
  !results printing loop
  Do N=0, NC
    YEX=FX(X(N+1))
    EF=DABS(YEX-Y(N+1))
    WRITE(*,100)  X(N + 1), Y(N + 1), YEX, EF, NL(N+1)
  End do
  PRINT *,'---------------------------------------------------------'  

10  format(' Starting X...: ')
20  format(' Ending X.....: ')
30  format(' Value Y at X0: ')
40  format(' Initial step : ')
50  format(' Desired error: ')
100 format(F6.1,F12.6,F12.6,F14.8,I6)

END

!Y' = F(X,Y)
real*8 Function F(XX,YY)
real*8 XX,YY
  F = XX * XX * YY
End

!Exact solution FX(XX)
real*8 Function FX(XX)
real*8 XX
  FX = DEXP(XX**3/3)
End

Subroutine Euler_Romberg(NC, H,ER,X0,Y0, X,Y, NL)
!***************************************************************
!* Solve an Ordinary Differential Equation Y'=F(X,Y) using the *
!* Euler-Romberg Method.                                       *
!* ----------------------------------------------------------- *
!* Inputs:                                                     *
!*          NC: number of points to calculate                  *
!*           H: abscissa integration step                      *
!*          ER: desired precision                              *
!*          X0: starting abscissa                              *
!*          Y0: initial value of Y at X=X0                     *
!*           F: external user defined function Y'              *
!* Outputs:                                                    *
!*           X: Table of NC abscissas                          *
!*           Y: Table of NC solution ordinates                 *
!*          NL: Table of number of steps per abscissa          *
!*                                                             *
!* Note:    the ending abscissa X1 is given by X0 + NC*H.      *
!***************************************************************                                                  
real*8 H,ER,X0,Y0,X(*),Y(*)
integer NL(*)
real*8 ET,XC,YC
real*8 T(20), F

  !Initial conditions
  X(0) = X0; Y(0) = Y0
  LA=10 !Maximum number of subdivisions
  !main integration loop
  DO N = 0, NC
    XC = X(N); YC = Y(N)
    T(1) = Y(N) + H * F(XC,YC)
    L = 1; LM = 2; ET=1.d0
    DO WHILE (L < LA.AND.ET >= ER)
      XC = X(N); YC = Y(N)
      DO J = 1, LM
        XC = XC + H / LM
        YC = YC + H / LM * F(XC,YC)
      END DO
      T(L + 1) = YC; M = 1; K = L; MM = 2; ET = 1.d0
      IF (K > 1) THEN
        DO WHILE (ET >= ER.AND.K > 1)
          T(K) = (MM * T(K + 1) - T(K)) / (MM - 1)
          ET = DABS(T(K) - T(K - 1))
          M = M + 1; K = K - 1; MM = MM * 2
        END DO
      END IF
      IF (K==1) THEN
        L = L + 1; LM = LM * 2
      END IF
    END DO
    X(N + 1) = X(N) + H; Y(N + 1) = T(K)
    NL(N+1)=L
  END DO
  Return  
End

!end of file eulromb1.f90