# Autocatalator.ode a' = mu*(kappa+c) - a*b^2 - a b' = (a*b^2 + a - b)/sigma c' = (b-c)/delta param mu=0.1, kappa=65, delta=2e-2, sigma=5e-3 # The variables range over several orders of magnitude, so # it's convenient to plot log(a), log(b) and log(c) aux la = log(a) aux lb = log(b) aux lc = log(c) # In order to avoid problems with the logs, # start from a point other than (0,0,0). a(0) = 1 b(0) = 1 c(0) = 1 # This system is stiff, so we need # an appropriate integrator. @ METHOD=stiff # The time scale of the oscillations is really fast, # and the spikes are really sharp and high, so we need # to adjust both the integration step size and # the maximum variable value allowed. @ DT=1e-4, BOUNDS=1e4 @ MAXSTOR=1000000 done