' PROGRAM WRITTEN BY INGEMAR BJERLE MARCH 2002
' SOLVING OF SECCOND ORDER LINEAR DIFFERENTIAL EQUATION
' INITIAL BOUNDARY VALUE SITUATION IS HANDLED

    [START]
      GOSUB [ini]
            INPUT "NUMBER OF STEPS (>100) ";N

            print "     XIN = ";XIN
            print "    XEND = ";XEND
            print "     YIN = ";YIN
            print "     UIN = ";UIN
            print "    XEND = ";XEND1
        PRINT""
        PRINT"    X         Y        Y1(analyt)"
        NN=10
      FOR K=0 TO NN
        DX=(XEND1-XIN)/N
        X=XIN:  Y=YIN:  U=UIN
          FOR I=1 TO N
            GOSUB [RUNGE]
          NEXT I
        Y1= X*X+.4*X+.08' ANALYTICAL SOLUTION
        PRINT USING("###.###",X);"   ";
        PRINT USING("#####.####",Y);"       ";Y1
        XEND1=XEND1+XEND/NN
      NEXT K
        PRINT""
        PRINT" FOLLOWS CLOSELY UP TO X=3, WITH 200 STEPS."
        PRINT" THE EQ STABLE AT YIN=.08. A SMALL DEVIATION FROM"
        PRINT" THE ANALYTICAL CURVE MAGNIFIES THE DEV EXPONENTIALLY."
        PRINT" RUNGE CAN NOT HANDLE THE SITUATION, EVEN WITH MANY STEPS"
        PRINT""
        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:  DY/DX - 5*(Y + X*X): FIRST ORDER ONLY USE A BELOW
        REM  WRITTEN AS TWO FIRST ORDER:
        REM  BOUNDARY CONDITIONS:  SEE [ini]
        'A=DY   B=DU  D=DZ
       '***************************CHOOSE THE SAME VARIABLE NAMES AS IN RUNGE
         A=5*(Y-X*X)
        'B=
        'D=
       '***************************INSERT THE BOUNDARY CONDITIONS IN MAIN PROGR
        RETURN

     [ini]     'SUBROUTINE WITH INITIAL DATES
          XIN=0
         XEND=5
          YIN=.08
          UIN=0
        XEND1=0
        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