'Program HENONLE.BAS calculates Lyapunov Exponent for Henon Map
'with B = 0.3 for values of C between 0 and 2.
'(c) 1997 by J. C. Sprott
SCREEN 12
DEFDBL A-Z
d0 = 1E-12
LINE (0, 0)-(639, 479), , B
LINE (0, 240)-(639, 240)
FOR i% = 1 TO 639
    c = i% / 320
    x1 = 0: y1 = 0
    x2 = 0: y2 = d0
    lsum = 0
    FOR j& = 1 TO 32100
        x1new = 1 - c * x1 * x1 + .3 * y1
        y1 = x1
        x1 = x1new
        x2new = 1 - c * x2 * x2 + .3 * y2
        y2 = x2
        x2 = x2new
        d1 = SQR((x2 - x1) ^ 2 + (y2 - y1) ^ 2)
        x2 = x1 + (d0 / d1) * (x2 - x1)
        y2 = y1 + (d0 / d1) * (y2 - y1)
        IF j& > 100 THEN lsum = lsum + LOG(d1 / d0)
        IF ABS(x1) > 1000 THEN END  'unbounded orbit
    NEXT j&
    le = lsum / (j& - 101)
    IF i% = 1 THEN PSET (i%, 240 - 480 * le) ELSE LINE -(i%, 240 - 480 * le)
    IF i% = 448 THEN LOCATE 2, 3: PRINT "LE ="; CSNG(le); "at c = 1.4"
NEXT i%
BEEP