CLS

DIM Y(5)
DIM Y0(5)
DIM W(5)
DIM Z(5)
DIM S(5)
DIM Q1(5)

REM ==RUNGE-KUTTA METHOD FOR A FIRST ORDER EQUATION==
REM READ DATA
GOSUB 150
X0 = A
H = H0
E5 = E1 / 100
50 :
PRINT
PRINT "> X = "; X0
FOR K = 1 TO M
PRINT "Y("; K; ") = "; Y0(K)
NEXT K
PRINT
GOSUB 650
REM ==CHECK FOR END INTERVAL==
IF (X0 - H) < B THEN 50
PRINT ">END OF EXECUTION"
END
150 :
PRINT " RUNGE-KUTTA METHOD FOR A"
PRINT " FIRST ORDER EQUATION"
PRINT
PRINT
PRINT ">NUMBER OF EQUATIONS"
PRINT
INPUT "M = "; M
PRINT
PRINT ">ENTER LIMITS OF INTEGRATION"
PRINT
INPUT " A = "; A
PRINT
INPUT " B = "; B
PRINT
PRINT ">ENTER INITIAL DATA"
PRINT
FOR K = 1 TO M
PRINT " Y"; K; "(A)= ";
INPUT Y0(K)
PRINT
NEXT K
PRINT ">INITIAL STEPSIZE "
PRINT
INPUT " H = "; H0
PRINT
PRINT ">LOCAL ERROR TOLERANCE "
PRINT
INPUT " E = "; E1
PRINT
PRINT ">STEPSIZE CONTROL? (Y/N) ";
INPUT C$
PRINT
PRINT
PRINT "SOLUTION"
PRINT "********"
PRINT
RETURN
650 :
T2 = H
T = T2
GOSUB 1100
FOR K = 1 TO M
Q1(K) = W(K)
NEXT K
710 :
T = T2 / 2
GOSUB 1100
FOR K = 1 TO M
Z(K) = W(K)
S(K) = Z(K)
NEXT K
X = X0 + T
GOSUB 1140
I8 = 0
FOR K = 1 TO M
Q3(K) = (W(K) - Q1(K)) / 15
D = ABS(Q3(K)) / T2
IF D >= E1 THEN 960
IF D < E5 THEN 860
I8 = 1
860 :
NEXT K
REM ==CONVERGENCE ATTAINED==
FOR K = 1 TO M
Y0(K) = W(K) + Q3(K)
NEXT K
H = T2
X0 = X0 + H
REM ==CHECK FOR STEP-SIZE CONTROL== IF C$ = "N" THEN 950
IF I8 = 1 THEN 950
H = 2 * H
950 :
RETURN
REM ==CONVERGENCE NOT ATTAINED==
960 :
IF C$ = "N" THEN 1020
T2 = T2 / 2
FOR K = 1 TO M
Q1(K) = Z(K)
NEXT K
GOTO 710
1020 :
PRINT
PRINT "**ERROR** CONVERGENCE "
PRINT "NOT ATTAINED WITHIN"
PRINT "REQUIRED TOLERANCE"
PRINT "FOR GIVEN STEP-SIZE"
PRINT "H="; H; " AT X= "; X0
END
REM ==RUNGE-KUTTA MODULE==
REM ==INITIALIZE VARIABLES==
1100 :
X = X0
FOR K = 1 TO M
S(K) = Y0(K)
NEXT K
1140 :
FOR K = 1 TO M
Y(K) = S(K)
W(K) = S(K)
NEXT K
REM ==COMPUTE W VALUE==
K = 0
1190 :
K = K + 1
GOSUB 2000
D = T * F
W(K) = W(K) + D / 6
Y(K) = S(K) + D / 2
IF K < M THEN 1190
X = X + T / 2
K = 0
1270 :
K = K + 1
GOSUB 2000
D = T * F
W(K) = W(K) + D / 3
Y(K) = S(K) + D / 2
IF K < M THEN 1270
K = 0
1340 :
K = K + 1
GOSUB 2000
D = T * F
Y(K) = S(K) + D
W(K) = W(K) + D / 3
IF K < M THEN 1340
X = X + T / 2
K = 0
1420 :
K = K + 1
GOSUB 2000
D = T * F
W(K) = W(K) + D / 6
IF K < M THEN 1420
RETURN

2000 :
IF K = 1 THEN 2010
IF K = 2 THEN 2030
IF K = 3 THEN 2050
IF K = 4 THEN 2070
IF K = 5 THEN 2090
2010 :
F = Y(2)
RETURN
2030 :
F = -.2 * Y(2)
RETURN
2050 :
F = Y(4)
RETURN
2070 :
F = -9.81 - .2 * Y(4)
RETURN
2090 :
F = 0
RETURN