'Program ABSCHAOS.BAS makes bifurcation diagram for chaotic flow with abs term
'(c) 1998 by J. C. Sprott <http://sprott.physics.wisc.edu/>
'Compile with PowerBASIC for DOS 3.5 <http://www.powerbasic.com/>
defext a-z
screen 12
amin=-.8
amax=-.5
xmin=-1
xmax=2
x=0: y=0: z=0
h=.002
tmax&=1000000
line(10,0)-(630,460),,b
locate 1,1: print trim$(str$(csng(xmax)));
locate 15,1: print"X";
locate 29,1: print trim$(str$(csng(xmin)));
locate 30,1: print trim$(str$(csng(amin)));
locate 30,41: print"A";
locate 30,81-len(trim$(str$(csng(amax)))): print trim$(str$(csng(amax)));
xticks%=6: yticks%=6
for i%=1 to xticks%-1
    xp%=10+620*i%/xticks%
    line(xp%,1)-(xp%,11)
    line(xp%,449)-(xp%,459)
next i%
for i%=1 to yticks%-1
    yp%=460-460*i%/yticks%
    line(11,yp%)-(21,yp%)
    line(619,yp%)-(629,yp%)
next i%
for i%=11 to 628
    a=amin+(amax-amin)*(i%-10)/620
    for t&=0 to tmax&
        call rk4(x,y,z,h)
        if abs(x)>1e3 then x=0: y=0: z=0: exit for
        if t&>tmax&/2 and x<xold and xold>xolder then
            call extremum (xolder, xold, x, dum, xe)
            j%=460-460*(xe-xmin)/(xmax-xmin)
            pset(i%,j%)
        end if
        xolder=xold
        xold=x
    next t&
    if inkey$=chr$(27) then end
next i%
while inkey$="": wend
end


def fnx(x,y,z)=y
def fny(x,y,z)=z
def fnz(x,y,z)=a*z-y+abs(x)-1

SUB rk4 (x, y, z, h)
    'Advances (x, y, z) with a fourth order Runge-Kutta method
    'User must provide the functions fnx, fny, fnz, and a step size h
    k1x = h * fnx(x, y, z)
    k1y = h * fny(x, y, z)
    k1z = h * fnz(x, y, z)
    k2x = h * fnx(x + .5 * k1x, y + .5 * k1y, z + .5 * k1z)
    k2y = h * fny(x + .5 * k1x, y + .5 * k1y, z + .5 * k1z)
    k2z = h * fnz(x + .5 * k1x, y + .5 * k1y, z + .5 * k1z)
    k3x = h * fnx(x + .5 * k2x, y + .5 * k2y, z + .5 * k2z)
    k3y = h * fny(x + .5 * k2x, y + .5 * k2y, z + .5 * k2z)
    k3z = h * fnz(x + .5 * k2x, y + .5 * k2y, z + .5 * k2z)
    k4x = h * fnx(x + k3x, y + k3y, z + k3z)
    k4y = h * fny(x + k3x, y + k3y, z + k3z)
    k4z = h * fnz(x + k3x, y + k3y, z + k3z)
    x = x + (k1x + 2 * (k2x + k3x) + k4x) / 6
    y = y + (k1y + 2 * (k2y + k3y) + k4y) / 6
    z = z + (k1z + 2 * (k2z + k3z) + k4z) / 6
END SUB

SUB extremum (yl, yo, yh, xe, ye)
    'Given three values yl @ x=-1, yo @ x=0 and yh @ x=1
    'Calculate the extremum ye and the value xe at which ye occurs
    'Fits the data to a parabola y = yo + bx + cx^2
    b = (yh - yl) / 2
    c = (yh + yl) / 2 - yo
    IF c = 0 THEN        'Second derivative is zero (no extremum)
    xe = 0
    ye = yo
    ELSE
    xe = -b / (2 * c)
    ye = yo - b * b / (4 * c)
    END IF
END SUB