      SUBROUTINE INITAL
C...
C...  EINSTEIN EQUATION FOR THE SCALE FACTOR OF THE UNIVERSE
C...
C...  THE FOLLOWING ORDINARY DIFFERENTIAL EQUATION (ODE) WAS PROPOSED
C...  BY EINSTEIN FOR THE SCALE FACTOR, A(T), OF THE EXPANSION OF THE
C...  UNIVERSE (1)
C...
C...      2    2
C...     D A/DT  = 8*PI*G*RHO*A/3                               (1)
C...
C...  WHERE
C...
C...     A(T)         SCALE FACTOR
C...
C...     T            TIME
C...
C...     G            NEWTON*S GRAVITATIONAL CONSTANT
C...
C...     RHO          VACUUM DENSITY
C...
C...  (1)  LAKE, GEORGE, WINDOWS ON A NEW COSMOLOGY, SCIENCE, V. 224,
C...  NO. 4650, PP 675-682, 18 MAY 85
C...
C...  EQUATION (1) IS INTEGRATED NUMERICALLY SUBJECT TO THE FOLLOWING
C...  INITIAL CONDITIONS FOR THE CURRENT EPOCH
C...
C...     A(0) = DA(0)/DT = 1                                 (2)(3)
C...
C...  AN ANALYTICAL SOLUTION TO EQUATIONS (1) TO (3) CAN EASILY BE
C...  DERIVED.  IF A SOLUTION OF THE FORM
C...
C...     A(T) = D*EXP(C*T)
C...
C...  IS ASSUMED, WHERE C = (8*PI*G*RHO/3)**(1/2), SUBSTUTION IN
C...  EQUATION (1) CONFIRMS THE SOLUTION.  THUS, THE GENERAL SOLUTION
C...  TO EQUATION (1) IS
C...
C...     A(T) = D1*EXP(C*T) + D2*EXP(-C*T)
C...
C...  WHERE D1 AND D2 ARE CONSTANTS TO BE DETERMINED FROM EQUATIONS (2)
C...  AND (3).  APPLYING THESE INITIAL CONDITIONS, WE HAVE
C...
C...     A(0) = D1 + D2 = 1
C...
C...     DA(0)/DT = C*D1 - C*D2 = 1
C...
C...  THUS, D2 = 1 - D1 AND
C...
C...     C*D1 - C*(1 - D1) = 1
C...
C...     2*C*D1 = 1 + C
C...
C...     D1 = (1 + C)/(2*C)
C...
C...  AND
C...
C...     D2 = 1 - (1 + C)/(2*C) = (C - 1)/(2*C)
C...
C...  THE SOLUTION IS THEREFORE
C...
C...     A(T) = (1 + C)/(2*C)*EXP(C*T) + (C - 1)/(2*C)*EXP(-C*T)   (4)
C...
C...  THE SOLUTION IS VERIFIED AS
C...
C...      2    2
C...     D A/DT  =  (C/2)*(1 + C)*EXP(C*T) - (C/2)*(1 - C)*EXP(-C*T)
C...
C...   -(C**2)*A = -(C/2)*(1 + C)*EXP(C*T) + (C/2)*(1 - C)*EXP(-C*T)
C...
C...     A(0) = (1 + C)/(2*C) - (1 - C)/(2*C) = 1
C...
C...     DA(0)/DT = (1/2)*(1 + C) + (1/2)*(1 - C) = 1
C...
C...  EQUATION (4) CAN THEN BE USED TO CHECK THE NUMERICAL SOLUTION.
C...
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/        A1,        A2
     2      /F/       A1T,       A2T
     3      /C/         C,        IP
     4      /H/         H
C...
C...  SET THE MODEL PARAMETERS
      C=1.
      C=SQRT(C)
C...

C...  INITIAL CONDITIONS
      A1=1.0
      A2=1.0
C...
C...  COMPUTE THE INITIAL DERIVATIVES
      CALL DERV
      IP=0
      RETURN
      END
      SUBROUTINE DERV
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/        A1,        A2
     2      /F/       A1T,       A2T
     3      /C/         C,        IP
     4      /H/         H
C...
C...  EQUATION (1) IS PROGRAMMED AS TWO FIRST-ORDER ODES
      A1T=A2
      A2T=C*A1
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/        A1,        A2
     2      /F/       A1T,       A2T
     3      /C/         C,        IP
     4      /H/         H
C...
C...  DIMENSION THE ARRAYS FOR PLOTTING
      DIMENSION HP(5), EP(5)
C...
C...  PRINT A HEADING FOR THE OUTPUT
      IP=IP+1
      IF(IP.EQ.1)THEN
         WRITE(NO,2)
2        FORMAT(1H ,//,' NORUN',9X,'T',9X,'H',8X,'A(T)',3X,
     1                    'ANAL A(T)',3X,'LOG ERROR')
      END IF
C...
C...  COMPUTE THE ANALYTICAL SOLUTION
      AE=(1.0+C)/(2.0*C)*EXP(C*T)-(1.0-C)/(2.0*C)*EXP(-C*T)
      ERR=(A1-AE)/AE*100.0
C...
C...  PRINT THE NUMERICAL AND ANALYTICAL SOLUTIONS
      WRITE(NO,1)NORUN,T,H,A1,AE,ERR
1     FORMAT(I6,2F10.4,3F12.4)
C...
C...  STORE THE ERROR IN THE NUMERICAL SOLUTION AND THE INTEGRATION
C...  STEP SIZE FOR SUBSEQUENT PLOTTING TO DETERMINE THE ORDER OF THE
C...  NUMERICAL METHOD (AT THE SIXTH POINT ALONG THE SOLUTION, WHICH
C...  IS SELECTED ARBITRARILY)
      IF(IP.EQ.6)THEN
         HP(NORUN)=ALOG10(ABS(H))
         EP(NORUN)=ALOG10(ABS(A1-AE))
      END IF
C...
C...  AT THE END OF THE FIFTH RUN, PRINT AND PLOT RESULTS FOR THE FIVE
C...  RUNS
      IF((NORUN.EQ.5).AND.(IP.EQ.11))THEN
C...
C...     PLOT THE ERROR VS INTEGRATION STEP
         CALL SPLOTS(1,NORUN,HP,EP)
         WRITE(NO,11)
11       FORMAT(1H ,//,' LOG(ERROR) VS LOG(H)')
C...
C...     COMPUTE AND PRINT THE ESTIMATED ORDER OF THE NUMERICAL ALGORITHM
C...     USING THE FIRST AND FIFTH POINTS
         ORDER=(EP(5)-EP(1))/(HP(5)-HP(1))
         WRITE(NO,12)ORDER
12       FORMAT(1H1,//,' ORDER OF METHOD = ',F5.3,' PTS 1 AND 5')
C...
C...     COMPUTE AND PRINT THE ESTIMATED ORDER OF THE NUMERICAL ALGORITHM
C...     USING THE THIRD AND FIFTH POINTS
         ORDER=(EP(5)-EP(3))/(HP(5)-HP(3))
         WRITE(NO,13)ORDER
13       FORMAT(1H1,//,' ORDER OF METHOD = ',F5.3,' PTS 3 AND 5')
      END IF
      RETURN
      END

