This is the BASIC code for solving a nth order IVP using the RK 4 method:
CLS:KEY OFF: "Program RK-n for nth order IVPs"
DIM A(4),G(4,4),W(4),X(4):W(1) = W(4) = 1/6:W(2)= W(3)=1/3
DEF FNF(T,X1,X2,X3,X4)= -4*SIN(T)- X1 -X2 - X3 + 4*X4 + X4*(X4- X3 + X2 - X1)
DEF FNQ(Z) = Z
INPUT"t0,n,h,kmax";T,N,H,KM
FOR I=1 TO N-1: PRINT "INPUT D";I;"(X)":INPUT X(N-I):NEXT I
FOR K =1 TO KM
PRINT T;X(N);EXP(T) + SIN(T)
FOR I=1 TO 4: FOR J = 1 TO N
P= .5:IF I =1 THEN P= 0: IF I = 4 THEN P = 1
IF J > 1 THEN G(I,J) = FNQ(X(J-1))
X(J)=X(J) + H*W(I)*G(I,J)
A(J)=X(J)+ P*H*G(I,J)
NEXT J
G(I,1)=FNF(T+P*H,A(1),A(2),A(3),A(4))
NEXT I:T = T + H:NEXT K: END
The problem is how is X(1) or A(1) calculated i.e. what happens when J = 1 ? This part :
IF J > 1 THEN G(I,J) = FNQ(X(J-1))
X(J)=X(J) + H*W(I)*G(I,J)
A(J)=X(J)+ P*H*G(I,J)
only gets run when J > 1 and this part:
G(I,1)=FNF(T+P*H,A(1),A(2),A(3),A(4))
requires A(1) to be calculated.
X(J)=...
andA(J)=...
lines are executed for every value ofJ
.THEN G(I,J) = FNQ(X(J-1))
is run only when J > 10
? (It's been a while since I used BASIC, and this might even depend on the implementation of it..... )