      PROGRAM HEUN

      INTEGER I, N
      REAL    H, X0, Y0, XFIM, X, Y
      REAL    K1, K2

      OPEN(UNIT = 1, FILE = 'heun.dat', STATUS = 'REPLACE')

      X0 = 0.0;
      Y0 = 1.0;

      H = 0.1;

      XFIM = 3.0;

      N = (XFIM - X0)/H;

      X = X0;
      Y = Y0;

      WRITE(1, *) X, Y

      DO 1 I = 0, N
           K1 = F(X, Y)

           K2 = F(X + H, Y + H * K1)

           Y  = Y + (H/2.0) * (K1 + K2)

           X = X + H

           WRITE(1, *) X, Y
1     CONTINUE

      CLOSE(1)
      STOP
      END

      REAL FUNCTION F(X, Y)
      F = Y
      RETURN
      END
