# stdliap.ode
# the liapunov exponent of the standard map
f(x)=mod(x+b+a*sin(2*pi*x),1)
fp(x)=1+2*pi*cos(2*pi*x)
init x=.1,a=0,z=0
par b=.75
x'=f(x)
a'=a
z'=z+heav(t-100)*ln(abs(fp(x))+1e-8)
aux liap=z/max(t-100,1)
@ xlo=2,xhi=4,ylo=-2,yhi=1
@ total=2000,trans=2000
@ meth=disc,bound=1000000
@ xp=a,yp=liap
done