;
;  program inte1.pro in IDL language 
;  (there is a free IDL called GDL you can download)
;
; solves numerically dx/dt = x(1-x) from t=0 to t=9
; to illustrate a second-order accurate integration scheme
;
; analytical solution to    dx/dt = x(1-x) 
; is, as you should be able to derive, 
; x/(1-x) = x0/(1-x0) * exp(t), where x0 is the starting x(t=0) 
; this can be re-written as 
; x(t) =  1 -  1/[1 + x0/(1-x0)*exp(t)]
; (check that x(t=0) indeed equals x0).
;

print,'      #steps           dt            T            x(T)          rel. error'
for j=1,12 do begin
t = 0d0  		; starting condition. 1d0 = 1 x 10^0 in double prec.
x = 0.1d0		; our startng condition
x_9 = 1d0-1d0/(1d0 +x/(1d0-x)*exp(9d0)) ; nearly dp-accurate result at t=9
thats_it = 0    ; a flag that will be 1 when integration reaches t=9    

dt0 = 10d0 ^ (-j/2d0)		; different timesteps will be studied
dt = dt0

; time integration loop
for i = 1L, 99000000L do begin      ; L means long integer (4B)
	if( (t+dt) GT 9d0) then begin   ; shorten dt if timestep goes to t>9 
		dt =  9d0 - t
		thats_it = 1  				; flag the end of the time loop
	endif

	t = t+dt
	f0 = x*(1d0-x)				; r.h.s. at the beg. of time step 
	x1 = x + f0*dt				; estimated x at the end of step	
	f1 = x1*(1d0-x1)			; estimated r.h.s. at the end of step
	x = x +dt*0.5d0*(f0+f1)     ; should be a 2nd order accurate new x

if(thats_it EQ 1) then begin
	print,i,dt0,t,x, (x/x_9 - 1.d0)  ; last column = relative accuracy
	break ; goto, here
	endif
endfor

here:

endfor
end

;IDL> .r inte1
;% Compiled module: $MAIN$.
;      #steps           dt            T            x(T)        rel. error ~ dt^2
;          29      0.31622777       9.0000000      0.99874769  -0.00014301319
;          91      0.10000000       9.0000000      0.99887822  -1.2332957e-05
;         285     0.031622777       9.0000000      0.99888936  -1.1826241e-06
;         901     0.010000000       9.0000000      0.99889043  -1.1714823e-07
;        2847    0.0031622777       9.0000000      0.99889053  -1.1673742e-08
;        9000    0.0010000000       9.0000000      0.99889054  -1.1661306e-09
;       28461   0.00031622777       9.0000000      0.99889054  -1.1656520e-10
;       90001   0.00010000000       9.0000000      0.99889054  -1.1645240e-11
;      284605   3.1622777e-05       9.0000000      0.99889054  -1.1773915e-12
;      900001   1.0000000e-05       9.0000000      0.99889054   5.9952043e-14 (*)
;     2846050   3.1622777e-06       9.0000000      0.99889054   2.8377301e-13 (*)
;     9000000   1.0000000e-06       9.0000000      0.99889054  -1.1324275e-14 (*)
;
; (*)- error governed by arithmetic accuracy and roundoff, not the integrator
;