q1'=(1-a1)*q1-a1*q2/2+a1*(M-c)/2
q2'=(1-a2)*q2-a2*q1/2+a2*(M-c)/2

init q1=1, q2=3
par  a1=0.1, a2=0.2, M=10, c=5
@ maxstor=100000, meth=discrete
done