function cnew = Reaction(cvec,nvec,dt,nSteps)
%----------------------------------------------------------------------------------- 
% Time marching the reaction step using the Crank-Nicholson rule
% @ Sohan Dharmaraja MIT JUN2007

coeff = 1;

for i=1:nSteps
    cnew = cvec + coeff*dt*(-nvec.*cvec) + (coeff^2*(dt^2)/2) * (nvec.^2).^cvec;
    cvec = cnew;
end