# the famous Lorenz equation set up for 3d view
init x=-7.5  y=-3.6  z=30
par r=27  s=10  b=2.66666  
# evaluate xdot etc
xd=s*(-x+y)
yd=r*x-y-x*z
zd=-b*z+x*y
# perturb x etc
xp=x+eps*vx
yp=y+eps*vy
zp=z+eps*vz
# evaluate perturbed xdot etc
xpd=s*(-xp+yp)
ypd=r*xp-yp-xp*zp
zpd=-b*zp+xp*yp
# approximate A(t) v
dx=(xpd-xd)/eps
dy=(ypd-yd)/eps
dz=(zpd-zd)/eps
# v^T A v
a=vx*dx+vy*dy+vz*dz
# odes for main system
x'=xd
y'=yd
z'=zd
# odes for vector
vx'=dx-a*vx
vy'=dy-a*vy
vz'=dz-a*vz
# compute avarage
aa'=a
aux abar=aa/max(1,t)
# check length
aux ll=vx^2+vy^2+vz^2
init vx=1
par eps=.001
@ dt=1, total=500,meth=qualrk,atol=1e-5,tol=1e-8
@ maxstor=20000,bound=1000000
done