' THIS PROGRAM MAKES THE SAME JOB AS DIFFEQ7.BAS
' PROGRMA WRITTEN BY INGEMAR BJERLE MARCH 2002
' SOLVING OF SECCOND ORDER  LINEAR DIFFERENTIAL EQUATION
' DIFF. EQ.:D2Y/DX2+Y/4-8=0
' TWO POINT BOUNDARY VALUE SITUATION IS HANDLED

      [START]
          INPUT "   STEP IN RUNGE (20-50) = ";N

          'SEARCHING PROCEDURE TO FIND DY/DX AT THE INLET
          'STARTING POINT REACTOR END
           XIN=0
          XEND=10
           YIN=0
           YEND=0
           UIN=0
           UEND=12*(-1)'NARROW GUESS NESSECARY
           DU=.01  'STARTING STEP SIZE IN SEARCH ROUTINE
            PRINT""
         PRINT"    X               Y             DY/DX"
          DX=(XIN-XEND)/N

          while DU>10^(-10)
          X=XEND:  Y=YEND:  U=UEND
           QQ=UEND
          FOR I=1 TO N
             GOSUB [RUNGE]
          NEXT I
          IF ABS(Y)>10^(-2) THEN
            UEND=UEND+DU
          ELSE
            UEND=UEND-DU:  DU=0.5*DU:  UEND=UEND+DU
          END IF
         WEND
          PRINT "    ";0;"             ";0;"            "; U;"  DY/DX AT INLET   "

           XIN=0
          XEND=10
           YIN=0
           YEND=0
           UIN=U  'HERE THE RESULT FROM THE SEARCHING PROCEDURE IS FED

        'THE PROBLEM FROM NOW ON IS AN INITIAL VALUE PROBLEM
        'AND THE SOLUTION IS STRAIGHT FORWARD
         NN=10
         DX=(XEND-XIN)/N/NN
         FOR KA=1 TO NN
          X=XIN:  Y=YIN:  U=UIN: Z=ZIN
          FOR I=1 TO N
             GOSUB [RUNGE]
          NEXT I
          PRINT"   ";X;"       ";Y;"        ";U
          XIN=XIN+XEND/NN
           YIN=Y
           UIN=U
           NEXT KA
           PRINT""
        'THE SOLUTION IS SYMMETRIC AT X=5. X=5 COULD HAVE BEEN CHOSEN AS
        ' STARTING POINT WITH KNOWN Y-VALUE AT X=5. AND WE WOULD HAVE HAD 
        AN INTIAL VALUE PROBLEM

        CONFIRM" DO YOU WANT TO QUIT Y/N";answer$
        IF answer$="yes" THEN
        GOTO [END]
        ELSE
        GOTO [START]
        END IF
  [END] PRINT " ENTER TO END"
        INPUT R$ 
        PRINT " *** END ***"
        END

     [EQ]       '  SUBROUTINE EQUATIONS
        REM  DIFF EQ:  D2Y/DX2+Y/4-8=0
        REM  WRITTEN AS TWO FIRST ORDER:  U=DY/DX:     DU/DX=8-Y/4
        REM  BOUNDARY CONDITIONS:  XIN=0:   XEND=10: Y(5)=71.9429: Y(10)=0
        REM  INTEGRATION STARTS AT END OF REACTOR. GIVES THE MOST STABLE SOLUTION
        'A=DY B=DU  D=DZ
       '*********************
        A=U
        B=8-Y/4
        'D=
       '*********************
     RETURN

     [INI]'NOT USED IN THIS CASE
     XIN=5
     XEND=10
     YIN=71.9429
     YEND=0
     UIN=0
     'UEND=
     RETURN

     [RUNGE]      ' SUBROUTINE RUNGE
        rem  VARIABLES IN SUBROUTINE: INDEPENTENT: X  DEPENDENT:  Y, U AND Z
        X=X: Y=Y: U=U: Z=Z
        GOSUB [EQ]
        K1=A*DX:    L1=B*DX:    F1=D*DX

        X=X+DX/2: Y=Y+K1/2: U=U+L1/2: Z=Z+F1/2
        GOSUB [EQ]
        K2=A*DX:    L2=B*DX:    F2=D*DX
        X=X-DX/2: Y=Y-K1/2: U=U-L1/2 :Z=Z-F1/2

        X=X+DX/2: Y=Y+K2/2: U=U+L2/2 :Z=Z+F2/2
        GOSUB [EQ]
        K3=A*DX:    L3=B*DX:    F3=D*DX
        X=X-DX/2: Y=Y-K2/2: U=U-L2/2 :Z=Z-F2/2

        X=X+DX:   Y=Y+K3:  U=U+L3    :Z=Z+F3
        GOSUB [EQ]
        K4=A*DX:    L4=B*DX:    F4=D*DX
        X=X-DX:   Y=Y-K3:  U=U-L3    :Z=Z-F3

        X=X+DX
        Y=Y+K1/6+K2/3+K3/3+K4/6:   U=U+L1/6+L2/3+L3/3+L4/6
        Z=Z+F1/6+F2/3+F3/3+F4/6

      RETURN