'Program FIG14.BAS shows an example of intermittency
'by J. C. Sprott <http://sprott.physics.wisc.edu/>
'Compile with the PowerBASIC Console Compiler <http://www.powerbasic.com/>
DEFEXT a-z      'do all calculations in extended (80-bit) precision

FUNCTION PBMAIN()
n&=3            'number of nonlinear odes
h=0.02          'stepsize of integration
tmax=4e3        'time to spend at each initial condition
s=100##         'scale factor on graph
pi=4#*ATN(1##)
CONSOLE NAME "Fig 14"
CONSOLE SET LOC 648,0
GRAPHIC WINDOW "Fig 14", 0, 0, 640, 480 TO hWin1&
GRAPHIC ATTACH hWin1&, 0
GRAPHIC WIDTH 1
GRAPHIC COLOR RGB(0,0,0), RGB(255,255,255)
GRAPHIC CLEAR
GRAPHIC ATTACH hWin1&, 0, REDRAW
GRAPHIC BOX(0,0)-(640,480)
GRAPHIC LINE(0,240)-(639,240)
DIM x(n&), dxdt(n&)
RANDOMIZE 13
x(1)=0.1*RND: x(2)=0.1*RND: x(3)=0.1*RND
PRINT x(1),x(2),x(3): WAITKEY$
FOR t=0 TO tmax
	FOR i&=1 TO 1/h
		CALL rk4(x(),dxdt(),n&,h,x())
		xp&=640*t/tmax
		yp&=240-240*x(1)/s
		GRAPHIC SET PIXEL(xp&,yp&)
	NEXT i&
	GRAPHIC REDRAW
NEXT t
GRAPHIC SAVE"fig14.bmp"
BEEP
WAITKEY$
END FUNCTION

SUB DERIVS (x(), dxdt(), n&)
	'Returns the time derivatives dxdt(i&) of x(i&)
	FOR i&=1 TO n&
		i1&=1+i& MOD n&
		dxdt(i&)=SIN(x(i1&))
	NEXT i&
END SUB

SUB RK4 (X(), DXDT(), N&, H, XNEW())
	'Fourth-order Runge-Kutta integrator
	DIM XT(N&), DXT(N&), DXM(N&)
	HH = H * .5##
	H6 = H / 6##
	CALL DERIVS(X(), DXDT(), N&)
	FOR I& = 1 TO N&
		XT(I&) = X(I&) + HH * DXDT(I&)
	NEXT I&
	CALL DERIVS(XT(), DXT(), N&)
	FOR I& = 1 TO N&
		XT(I&) = X(I&) + HH * DXT(I&)
	NEXT I&
	CALL DERIVS(XT(), DXM(), N&)
	FOR I& = 1 TO N&
		XT(I&) = X(I&) + H * DXM(I&)
		DXM(I&) = DXT(I&) + DXM(I&)
	NEXT I&
	CALL DERIVS(XT(), DXT(), N&)
	FOR I& = 1 TO N&
		XNEW(I&) = X(I&) + H6 * (DXDT(I&) + DXT(I&) + 2## * DXM(I&))
	NEXT I&
	ERASE DXM, DXT, XT
END SUB