100'  NAME--DE-SYS
110'
120'  DESCRIPTION--SOLVES THE INITIAL VALUE PROBLEM
130'     X' = F(T,X,Y)
140'     Y' = G(T,X,Y)
150'     X(T0) = X0
160'     Y(T0) = Y0
170'
180'  SOURCE--RUNGA KUTTA METHOD (2ND ORDER ACCURACY)
190'
200'  INSTRUCTIONS--THE INTEGRATION STEP IS H AND THE SOLUTION IS
210'  TABULATED ON THE INTERVAL   T0<= T <= B IN STEPS OF SIZE L.
220'
230'     390 DEF FNF(T,X,Y) = ......
240'     400 DEF FNF(T,X,Y) = ......
250'     410 DATA T0,X0,Y0,B
260'     420 DATA H,L
270'
280'  FOR MODIFICATIONS OF THE PRINT ROUTINE, FOR EXAMPLE TO
290'  PLOT THE SOLUTION RATHER THAN TABULATE ITS VALUES, THE 
300'  FOLLOWING LINES SHOULD BE MODIFIED:
310'     510 PRINT"VALUE OF T","VALUE OF X",VALUE OF Y"
320'     520 PRINT
330'     530 PRINT T,X,Y
340'     800 PRINT Q,X2,Y2
350'
360'
370'  *  *  *  *  *   MAIN PROGRAM  *  *  *  *  *  *  *  *  *  *  *  *
380'
390 DEF FNF(T,X,Y)=X-X*Y
400 DEF FNG(T,X,Y)=X*Y-Y
410DATA 0,.5,.5,4
420DATA .01,.1
430 READ T0,X,Y,B,H,L
440 LET T=T0
450 LET E=1/L
460 IF (B-T0)*H>0 THEN 480
470 LET H=-H
480 IF (B-T0)*L>0 THEN 500
490 LET L=-L
500 REM
510 PRINT "VALUE OF T","VALUE OF X","VALUE OF Y"
520 PRINT
530 PRINT T,X,Y
540 GOTO 600
550 LET P1=X+H*FNF(T,X,Y)
560 LET Q1=Y+H*FNG(T,X,Y)
570 LET X=(X+P1)/2+.5*H*FNF(T+H,P1,Q1)
580 LET Y=(Y+Q1)/2+.5*H*FNG(T+H,P1,Q1)
590 LET T=T+H
600 IF (T+H-B)*SGN(L)+1E-7 >=0 THEN 630
610 LET R=T+H
620 GOTO 640
630 LET R=B
640 LET A=INT(E*R)/E + L*(SGN(H)-1)/2
650 IF (A+L-R)*SGN(L)-1E-7 >= 0 THEN 670
660 LET A=A+L
670 IF (T-A)*SGN(L)+1E-7 >=0 THEN 700
680 LET Q=A
690 GOSUB 760
700 IF R=B THEN 720
710 GOTO 550
720 IF ABS(R-A)<.5*1E-6 THEN 750
730 LET Q=R
740 GOSUB 760
750 STOP
760 LET P2=X+(Q-T)*FNF(T,X,Y)
770 LET Q2=Y+(Q-T)*FNG(T,X,Y)
780 LET X2=(X+P2)/2+.5*(Q-T)*FNF(Q,P2,Q2)
790 LET Y2=(Y+Q2)/2+.5*(Q-T)*FNG(Q,P2,Q2)
800 PRINT Q,X2,Y2
810 RETURN
820 END