I am trying to numerically solve a coupled second order differential equation that describes the rate of a reversible equation between two substances A and B.

where k- and k+ are the rates of conversion during the reaction, D is the diffusion coefficient and v is the advection velocity.

The system is assumed to be in a steady state where dA/dt = dB/dt = 0.

I am trying to use the finite difference method to solve this, while my graph is convergent, it is not stable.

 /* Loop over points */
for (j=0; j<len; j++) {
  /* Centred space indecis for periodic boundary */
  jm = (j+len-1)%len;
  jp = (j+    1)%len;
  /* Finite difference evaulation of gradient */
  lslopeA = (A[jp] - A[jm]) / ( 2 * dx );
  d2Adx2 = (y[jp] + y[jm] - 2*y[j]) / (dxsq);
  lslopeB= (B[jp] - B[jm]) / ( 2 * dx );
  d2Bdx2 = (B[jp] + B[jm]] - 2*z[j]) / (dxsq);
  /* Forward time step */
  A_next[j] = A[j] - (lslopeA * ts * vel) + (d2Adx2 * D * ts) - ts * (kp * A[j] - km * B[j] - S[j]) ;
  B_next[j] = B[j] - (lslopeB * ts * vel) + (d2Bdx2 * D * ts) + ts * (kp * A[j] - km * B[j] + sigma[j] * B[j]) ;
}
/* Copy next time step solution to current one */
for (j=0; j<len; j++) {
  A[j] = A_next[j];
  B[j] = B_next[j];
}

I hope you can tell me where i'm going wrong, my values are all below the stability criteria set by the scheme.