'
' *****************************************************************
'
'                       FROM CALCULUS TO CHAOS
'       
'                    An Introduction to Dynamics
'
'                  
'                                by
'
'                           David Acheson
'
' ******************************************************************
'
'        The program which follows is part of the above book,
'            published in 1997 by Oxford University Press.
'                   ISBN 0-19-850077-7 (paperback)
'              
'
'                   Copyright ¸ David Acheson 1997
'
'
'
'
'                       Program 6:  NVARY.BAS
'
'
DEFDBL A-H, K-M, O-Z: DEFINT I-J, N
n = 3: OPTION BASE 1
DIM x(n), xc(n), f(n), c1(n), c2(n), c3(n), c4(n)
DIM xcold(n)

REM ****** Setting up graphics ******
 
 CLS : SCREEN 9: PAINT (1, 1), 9
  xm = 8: ym = 8
 VIEW (180, 17)-(595, 330), 0, 14
 WINDOW (-xm, -ym)-(xm, ym)
 LINE (-xm, 0)-(xm, 0), 11
 LINE (0, -ym)-(0, ym), 11
  LOCATE 13, 76: PRINT xm: LOCATE 1, 48: PRINT ym
  LOCATE 15, 1: PRINT "xi"
  LOCATE 23, 2: PRINT "Time"

REM ****** Step-by-step method ******
 
 LOCATE 13, 2: INPUT "w"; w

 k = .1#: b = .04#

 xc(1) = 0#: xc(2) = 0#: xc(3) = 0#
 h = .1#
   
 DO
  
   FOR i = 1 TO n: xcold(i) = xc(i): NEXT
        CLS
        LINE (-xm, 0)-(xm, 0), 11
        LINE (0, -ym)-(0, ym), 11
        t = 0#
    
     DO
       GOSUB Runge
       t = t + h
       PSET (xc(1), xc(2)), 14
       LOCATE 23, 6: PRINT INT(t)
     LOOP UNTIL INKEY$ = "q"
      LOCATE 16, 1
       FOR i = 1 TO n: PRINT "              ": NEXT
      LOCATE 16, 1
       FOR i = 1 TO n: PRINT CSNG(xc(i)): NEXT
              
   LOCATE 13, 2: INPUT "new w"; w
    IF w = 0# THEN
      FOR i = 1 TO n: xc(i) = xcold(i): NEXT
      LOCATE 13, 2: INPUT "new w"; w
    END IF

 LOOP

REM ****** Subroutines ******

Equations:
   
 f(1) = x(2)
 f(2) = -k * x(2) / w + (-x(1) - b * x(1) ^ 3# + COS(x(3))) / w ^ 2#
 f(3) = 1#
 RETURN

Runge:

   FOR i = 1 TO n: x(i) = xc(i): NEXT
   GOSUB Equations
   FOR i = 1 TO n: c1(i) = h * f(i): NEXT

   FOR i = 1 TO n: x(i) = xc(i) + c1(i) / 2#: NEXT
   GOSUB Equations
   FOR i = 1 TO n: c2(i) = h * f(i): NEXT

   FOR i = 1 TO n: x(i) = xc(i) + c2(i) / 2#: NEXT
   GOSUB Equations
   FOR i = 1 TO n: c3(i) = h * f(i): NEXT

   FOR i = 1 TO n: x(i) = xc(i) + c3(i): NEXT
   GOSUB Equations
   FOR i = 1 TO n: c4(i) = h * f(i): NEXT

 FOR i = 1 TO n
 xc(i) = xc(i) + (c1(i) + 2# * c2(i) + 2# * c3(i) + c4(i)) / 6#
 NEXT
   
 RETURN