CLS
DIM Y(5)
DIM Y0(5)
DIM W(5)
DIM Z(5)
DIM S(5)
DIM Q1(5)
REM ==RUNGE-KUTTA METH0D FOR A FIRST ORDER EQUATION==
REM READ DATA
GOSUB [ONFIZE]
X0 = A
H = H0
E5 = E1/100
[FIZE]
PRINT
PRINT "X = ";X0
FOR K = 1 TO M
PRINT "Y";K;" =";Y0(K)
NEXT K
PRINT
GOSUB [SIFIZE]
REM ==CHECK FOR END INTERVAL==
IF (X0-H) < B THEN [FIZE]
PRINT "END OF EXECUTION"
END
[ONFIZE]
PRINT "RUNGE-KUTTA METH0D 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 " ";YY
Y0(K) = YY
PRINT
NEXT K
PRINT "INITIAL STEPSIZE "
PRINT
INPUT "H =";H0
PRINT
PRINT "LOCAL ERROR TOLERANCE X AS IN 1E-X"
PRINT
INPUT "E = ";E1
E1=1/10^E1
PRINT
PRINT "STEPSIZE CONTROL?(Y/N)";
INPUT C$
PRINT
PRINT
PRINT "SOLUTION"
PRINT "********"
PRINT
RETURN
[SIFIZE]
T2 = H
T = T2
GOSUB [ONONZEZE]
FOR K = 1 TO M
Q1(K) = W(K)
NEXT K
[SEONZE]
T = T2 / 2
GOSUB [ONONZEZE]
FOR K = 1 TO M
Z(K) = W(K)
S(K) = Z(K)
NEXT K
X = X0 + T
GOSUB [ONONFOZE]
I8 = 0
FOR K = 1 TO M
Q3(K) = (W(K)-Q1(K))/15
D = ABS(Q3(K))/T2
IF D >= E1 THEN [NISIZE]
IF D < E5 THEN [EISIZE]
I8 = 1
[EISIZE]
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 [NIFIZE]
IF I8 = 1 THEN [NIFIZE]
H = 2 * H
[NIFIZE]
RETURN
REM ==CONVERGENCE NOT ATTAINED==
[NISIZE]
IF C$ ="N" THEN [ONZETWZE]
T2 = T2/2
FOR K = 1 TO M
Q1(K) = Z(K)
NEXT K
GOTO [SEONZE]
[ONZETWZE]
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 ==
[ONONZEZE]
X = X0
FOR K = 1 TO M
S(K) = Y0(K)
NEXT K
[ONONFOZE]
FOR K = 1 TO M
Y(K) = S(K)
W(K) = S(K)
NEXT K
REM ==COMPUTE W VALUE==
K = 0
[ONONNIZE]
K = K + 1
GOSUB [TWZEZEZE]
D=T *F
W(K) = W(K) + D/6
Y(K) = S(K) + D/2
IF K < M THEN [ONONNIZE]
X = X + T/2
K = 0
[ONTWSEZE]
K = K+ 1
GOSUB [TWZEZEZE]
D = T * F
W(K) = W(K) + D/3
Y(K) = S(K) + D/2
IF K < M THEN [ONTWSEZE]
K = 0
[ONTHFOZE]
K = K+1
GOSUB [TWZEZEZE]
D =T * F
Y(K) = S(K) + D
W(K) = W(K) + D/3
IF K < M THEN [ONTHFOZE]
X = X + T/2
K = 0
[ONFOTWZE]
K = K + 1
GOSUB [TWZEZEZE]
D = T * F
W(K) = W(K) + D/6
IF K < M THEN [ONFOTWZE]
RETURN
[TWZEZEZE]
IF K = 1 THEN [TWZEONZE]
IF K = 2 THEN [TWZETHZE]
IF K = 3 THEN [TWZEFIZE]
IF K = 4 THEN [TWZESEZE]
IF K = 5 THEN [TWZENIZE]
[TWZEONZE]
F = Y(2)
RETURN
[TWZETHZE]
F = -.2*Y(2)
RETURN
[TWZEFIZE]
F = Y(4)
RETURN
[TWZESEZE]
F = -9.81-.2*Y(4)
RETURN
[TWZENIZE]
F =  0