' PROGRAM WRITTEN BY INGEMAR BJERLE MARCH 2002
' SOLVING OF SECCOND ORDER NON LINEAR DIFFERENTIAL EQUATION
' THE EQUATION KNOWN AS THE DISPERSION MODEL BY CHEMICAL ENGINEERS
' TWO POINT BOUNDARY VALUE SITUATION IS HANDLED

      [START]
         INPUT" GAMMA(.1-5)= ";G   'DISPERSION PARAMETER RANGE: 0.1-5
         INPUT"    R(.5-50)= ";R   'REACTION PARAMETER   RANGE: 0.1-50
         INPUT"STEPS(1-100)= ";N   'STEPS IN RUNGE RANGE: 10-50
         PRINT " GAMMA   Y-IN         Y-END"

      for M=1 to 20          'GAMMA IS CHANGED IN STEPS OF 0.1
          YEND=EXP((-1)*R)   'ANALYTICAL END-VALUE FOR IDEAL TUBULAR REACTOR
          DY=0.2             'STARTING STEP SIZE IN SEARCH ROUTINE
          XIN=0
          XEND=1
          UEND=0
          DX=(XIN-XEND)/N
         while DY > 10^(-5)
          X=XEND:  Y=YEND:  U=UEND
          FOR I=1 TO N
             GOSUB [RUNGE]
          NEXT I
          IF ((Y-1)/G-U)>0 THEN
            YEND=YEND-DY:  DY=0.5*DY:  YEND=YEND+DY
          ELSE
            YEND=YEND+DY
          END IF
         WEND
          PRINT G;"  "; Y;"     "; YEND
          G=G+0.1
       NEXT M

        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:  G*D2Y/DX2 - DY/DX - R*Y^2
        REM  WRITTEN AS TWO FIRST ORDER:  U=DY/DX:     G*DU/DX-U-R*Y^2
        REM  BOUNDARY CONDITIONS:  XIN=0:   XEND=1:   UEND=0:  (YIN-1)/G-UIN)=0: TO FIND BY SEARCH
        REM  INTEGRATION STARTS AT END OF REACTOR. GIVES THE MOST STABLE SOLUTION
        REM  A=DY, B=DU,  D=DZ
       '*********************
        A=U
        B=(U+R*Y*Y)/G
        'D=
       '*********************
     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