        PROGRAM D15R2

C       Driver for routine RKDUMB

        PARAMETER(NVAR=2)

        DIMENSION VSTART(NVAR)

        COMMON /PATH/X(200),Y(10,200)

        OPEN(6,FILE='RK.OUT',FORM='FORMATTED')

        WRITE(*,*) ' X1=?  X2=?  NSTEP=?'

        READ(*,*) X1,X2,NSTEP

        VSTART(1)=1.

        VSTART(2)=0.

        CALL RKDUMB(VSTART,X1,X2,NSTEP)

        WRITE(*,'(/1X,T9,A,T17,A,T31,A/)') 'X','X1 VALUE','X2 VALUE'

        WRITE(6,'(/1X,T9,A,T17,A,T31,A/)') 'X','X1 VALUE','X2 VALUE'

        DO 11 I=1,NSTEP+1

                WRITE(*,'(1X,F10.4,2X,2F12.6)') X(i),Y(1,i),Y(2,i)

                WRITE(6,'(1X,F10.4,2X,2F12.6)') X(i),Y(1,i),Y(2,i)

11      CONTINUE

        END

C

        SUBROUTINE DERIVS(X,Y,DYDX)

        PARAMETER(NVAR=2)

        DIMENSION Y(NVAR),DYDX(NVAR)

        DYDX(1)=Y(2)

        DYDX(2)=-Y(2)*x-y(1)

        RETURN

        END

C

      SUBROUTINE RKDUMB(VSTART,X1,X2,NSTEP)

      PARAMETER (NMAX=10,NVAR=2)

      COMMON /PATH/ XX(200),Y(10,200)

      DIMENSION VSTART(NVAR),V(NMAX),DV(NMAX)

      DO 11 I=1,NVAR

        V(I)=VSTART(I)

        Y(I,1)=V(I)

11    CONTINUE

      XX(1)=X1

      X=X1

      H=(X2-X1)/NSTEP

      DO 13 K=1,NSTEP

        CALL DERIVS(X,V,DV)

        CALL RK4(V,DV,X,H,V)

        IF(X+H.EQ.X)PAUSE 'Stepsize not significant in RKDUMB.'

        X=X+H

        XX(K+1)=X

        DO 12 I=1,NVAR

          Y(I,K+1)=V(I)

12      CONTINUE

13    CONTINUE

      RETURN

      END

C

      SUBROUTINE RK4(Y,DYDX,X,H,YOUT)

      PARAMETER (NMAX=10,NVAR=2)

      DIMENSION Y(NVAR),DYDX(NVAR),YOUT(NVAR),YT(NMAX),DYT(NMAX),

     @DYM(NMAX)

      HH=H*0.5

      H6=H/6.

      XH=X+HH

      DO 11 I=1,NVAR

        YT(I)=Y(I)+HH*DYDX(I)

11    CONTINUE

      CALL DERIVS(XH,YT,DYT)

      DO 12 I=1,NVAR

        YT(I)=Y(I)+HH*DYT(I)

12    CONTINUE

      CALL DERIVS(XH,YT,DYM)

      DO 13 I=1,NVAR

        YT(I)=Y(I)+H*DYM(I)

        DYM(I)=DYT(I)+DYM(I)

13    CONTINUE

      CALL DERIVS(X+H,YT,DYT)

      DO 14 I=1,NVAR

        YOUT(I)=Y(I)+H6*(DYDX(I)+DYT(I)+2.*DYM(I))

14    CONTINUE

      RETURN

      END

