C     TITLE: DEMONSTRATION OF EULER'S METHOD FOR SYSTEMS.
C
C     THIS SOLVES THE INITIAL VALUE PROBLEM FOR THE FIRST ORDER
C     SYSTEM
C        Y'(X) = F(X,Y(X)) ,  X0 .LE. X .LE. B ,  Y(X0)=Y0.
C     THE FUNCTION F(X,Z) IS DEFINED BELOW, ALONG WITH THE TRUE
C     SOLUTION Y(X).  THE NUMBER OF THE PROBLEM TO BE SOLVED
C     IS SPECIFIED BY THE INPUT VARIABLE 'NUMDE', WHICH IS USED
C     IN THE SUBROUTINES 'FCN' AND 'TRUE'.  THE PROGRAM WILL
C     REQUEST THE PROBLEM PARAMETERS, ALONG WITH THE VALUES OF
C     'H' AND 'IPRINT'.  'H' IS THE STEPSIZE, AND 'IPRINT' IS
C     THE NUMBER OF STEPS BETWEEN EACH PRINTING OF THE SOLUTION.
C     USE H=0 AND NUMDE=0 TO STOP THE PROGRAM.
C
C     TO INCREASE THE ORDER OF THE SYSTEMS WHICH CAN BE
C     HANDLED, INCREASE THE SIZE OF 'MAXN' IN THE PARAMETER
C     STATEMENT GIVEN BELOW.
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (MAXN=3)
      DIMENSION Y0(MAXN), Y1(MAXN), Y(MAXN), YZERO(MAXN), F(MAXN) 
      PARAMETER (ZERO=0.0D0)
      COMMON/BLOCKF/NUMDE
C
C     INPUT PROBLEM PARAMETERS.
10    PRINT *, ' WHICH DIFFERENTIAL EQUATION?'
      PRINT *, ' GIVE ZERO TO STOP.'
      READ *, NUMDE
      IF(NUMDE .EQ. 0) STOP
      PRINT *, ' GIVE DOMAIN [X0,B] OF SOLUTION.'
      READ *, XZERO, XEND
      PRINT *, ' GIVE THE ORDER OF THE SYSTEM.  IT SHOULD BE'
      PRINT *, ' LESS THAN OR EQUAL TO MAXN=', MAXN
      READ *, N
      PRINT *, ' GIVE THE COMPONENTS OF Y0, ONE AT A TIME.'
      DO 15 J=1,N
        PRINT *, ' GIVE COMPONENT', J
15      READ *, YZERO(J)
C
20    PRINT *, ' GIVE STEPSIZE H AND PRINT PARAMETER IPRINT.'
      PRINT *, ' LET H=0 TO TRY ANOTHER DIFFERENTIAL EQUATION.'
      READ *, H,IPRINT
      IF(H .EQ. ZERO) GO TO 10
C
C     INITIALIZE.
      X0 = XZERO
      DO 25 J=1,N
25      Y0(J) = YZERO(J)
      PRINT 1000, NUMDE, XZERO, XEND, H, IPRINT
1000  FORMAT(//,' EQUATION',I2,5X,'XZERO =',1PE9.2,5X,'B =',E9.2,
     *       /,' STEPSIZE =',E10.3,5X,'PRINT PARAMETER =',I3)
      PRINT 1001, N
1001  FORMAT(' ORDER =',I1,//,' J',4X,'Y0(J)')
      PRINT 1002, (J,YZERO(J),J=1,N)
1002  FORMAT(I2,2X,1PE9.2)
      PRINT 1003
1003  FORMAT(/,5X,'X',7X,'J',6X,'APPROX Y(J)',7X,'ERROR')
C
C     BEGIN MAIN LOOP FOR COMPUTING SOLUTION OF DIFFERENTIAL EQUATION.
30    DO 40 K=1,IPRINT
        X1 = X0 + H
        IF(X1 .GT. XEND) GO TO 20
        CALL FCN(X0,Y0,F)
        DO 35 J=1,N
          Y1(J) = Y0(J) + H*F(J)
35        Y0(J) = Y1(J)
        X0 = X1
40      CONTINUE
C
C     CALCULATE ERROR AND PRINT RESULTS.
      CALL TRUE(X1,Y)
      DO 60 J=1,N
60      PRINT 1004, X1, J, Y1(J), Y(J) - Y1(J)
1004  FORMAT(F10.5,3X,I1,5X,1PE13.6,4X,E9.2)
      GO TO 30
      END

      SUBROUTINE FCN(X,Z,F) 
C     -------------- 
C
C     THIS DEFINES THE RIGHT SIDE OF THE
C     DIFFERENTIAL EQUATION, F(X,Z).
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z) 
      PARAMETER (TWO=2.0D0, THREE=3.0D0, FOUR=4.0D0, FIVE=5.0D0)
      DIMENSION Z(*), F(*)
      COMMON/BLOCKF/NUMDE
C
      GO TO (10,20), NUMDE
10    CS = COS(X)
      SN = SIN(X)
      F(1 )= Z(1) - TWO*Z(2) + FOUR*CS - TWO*SN
      F(2) = THREE*Z(1) - FOUR*Z(2) + FIVE*(CS - SN)
      RETURN
20    F(1) = Z(2)
      F(2) = Z(3)
      F(3) = -THREE*(Z(3) + Z(2)) - Z(1) - FOUR*SIN(X)
      RETURN            
      END

      SUBROUTINE TRUE(X,Y)  
C     ---------------
C
C     THIS GIVES THE TRUE SOLUTION OF
C     THE INITIAL VALUE PROBLEM, Y(X).
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z) 
      PARAMETER (TWO=2.0D0)
      DIMENSION Y(*)
      COMMON/BLOCKF/NUMDE
C
      GO TO(10,20), NUMDE
10    CS = COS(X)
      SN = SIN(X)
      Y(1) = SN + CS
      Y(2) = TWO*CS
      RETURN
20    SN = SIN(X)
      CS = COS(X)
      Y(1) = SN + CS
      Y(2) = CS - SN
      Y(3) = -SN - CS
      RETURN
      END