function xdot=reac1(t,x,a0,b0,k,r)
xdot=k*(a0-x/2)*(b0-x/2)-r*x^2;

