'
' *****************************************************************
'
'                       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 2:  2PHASE.BAS
'
'
REM ****** Setting up graphics ******

CLS : SCREEN 9: PAINT (1, 1), 9
xm = 4: ym = 4: tm = 25
VIEW (180, 17)-(595, 330), 0, 13
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

REM ****** Functions f(x,y) and g(x,y) ******

  DEF fnf (x, y) = y
  DEF fng (x, y) = -SIN(x)

REM ****** Direction field ******
  
p = 25
   FOR y = ym TO -ym STEP -2 * ym / p
     FOR x = -xm TO xm STEP 2 * xm / p
       x1 = fnf(x, y) / xm: y1 = fng(x, y) / ym
       s = p * SQR(x1 ^ 2 + y1 ^ 2)
       x2 = fnf(x, y) / s: y2 = fng(x, y) / s
       LINE (x, y)-(x + x2, y + y2), 9
       CIRCLE (x + x2, y + y2), .003 * xm, 9
     NEXT x
   NEXT y

REM ****** Step-by-step method ******

  DO
     t = 0: LOCATE 13, 1: INPUT "x0,y0"; x, y
     h = .01
    DO
       GOSUB ImpEuler
       t = t + h
       PSET (x, y), 14
    LOOP UNTIL INKEY$ = "q"
  LOOP

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

Euler:

    c1 = h * fnf(x, y)
    d1 = h * fng(x, y)
    x = x + c1
    y = y + d1
    RETURN

ImpEuler:

    c1 = h * fnf(x, y)
     d1 = h * fng(x, y)
    c2 = h * fnf(x + c1, y + d1)
     d2 = h * fng(x + c1, y + d1)
    x = x + (c1 + c2) / 2
     y = y + (d1 + d2) / 2
    RETURN