Wrong reaction order in an autocatalytic reaction

**Please complete the following tasks**

- [x] I have searched the documentation and help resources before reporting.
- [x] I have searched the [discussions](https://github.com/geem-lab/overreact/discussions) to avoid duplicates.
- [x] I have searched the [open](https://github.com/geem-lab/overreact/issues) and [closed](https://github.com/geem-lab/overreact/issues?q=is%3Aissue+is%3Aclosed) issues to avoid duplicates.

**Describe the bug**
Consider the autocatalytic reaction "2B -> B + C":

```python
import numpy as np

import overreact as rx

scheme = rx.parse_reactions("2 B -> B + C")

stoichiometry = np.asarray(scheme.A)
kinetic_orders = np.where(stoichiometry > 0, 0, -stoichiometry)
reaction_order = np.sum(kinetic_orders[:, 0])

print("reaction (should be 2nd order in B):", scheme.reactions[0])
print("computed net stoichiometry (A):", stoichiometry[:, 0])
print("reaction order:", reaction_order)

# actually computing the schema to show it's first order:
dydt = rx.get_dydt(scheme, k=np.array([1.0]))
rate_b1 = dydt(0.0, np.array([1.0, 0.0]))[1]
rate_b2 = dydt(0.0, np.array([2.0, 0.0]))[1]

print("rate at [B] = 1:", float(rate_b1))
print("rate at [B] = 2:", float(rate_b2))
```

Output:

```
reaction (should be 2nd order in B): 2 B -> B + C
computed net stoichiometry (A): [-1. 1.]
reaction order: 1.0
An NVIDIA GPU may be present on this machine, but a CUDA-enabled jaxlib is not installed. Falling back to cpu.
rate at [B] = 1: 1.0
rate at [B] = 2: 2.0
```

**To Reproduce**
Steps to reproduce the behavior:

1. Open a terminal.
2. Save the above in a script with a proper overreact environment.
3. Execute and see the inferred reaction order in B is 1 instead of 2.

**Expected behavior**
Should be 2nd order in B.

**Origin**

I think [_add_reaction()](https://github.com/geem-lab/overreact/blob/dec053bb108c117fef10eee446be06655384ecdc/overreact/core.py#L603-L639) when taken together with [_dydt()](https://github.com/geem-lab/overreact/blob/dec053bb108c117fef10eee446be06655384ecdc/overreact/simulate.py#L262-L265) is the faulty point.

```python
def _add_reaction(reactants, products, is_half_equilibrium, transition):
"""Local helper function with side-effects."""
# TODO(schneiderfelipe): what if reaction is defined, then redefined
# as equilibrium, or vice-versa?
if (
transition is not None
and (reactants, products, is_half_equilibrium, transition) in reactions
):
return

# found new reaction
reactions[(reactants, products, is_half_equilibrium, transition)] = None

A_vector = np.zeros(len(compounds))
for coefficient, reactant in reactants:
A_vector[compounds[reactant]] = -coefficient
B_vector = A_vector

if transition is not None:
B_vector = A_vector.copy()

# it's assumed that
# 1. there's a singe transition compound, and
# 2. its coefficient equals one
B_vector[compounds[transition[-1][-1]]] = 1

for coefficient, product in products:
A_vector[compounds[product]] += coefficient

if (
is_half_equilibrium
and (products, reactants, is_half_equilibrium, transition) in reactions
):
B_vector = np.zeros(len(compounds))

A.append(A_vector)
B.append(B_vector)
```

The parser sotres the reactants:

A_vector[B] = -2

Then it adds the products:

A_vector[B] += 1
A_vector[C] += 1

Which is fine for the net concentration changes, but when we get to [_dydt()](https://github.com/geem-lab/overreact/blob/dec053bb108c117fef10eee446be06655384ecdc/overreact/simulate.py#L262-L265), more specifically, in the line below:

```python
M = jnp.where(A > 0, 0, -A).T
```

we compute M based on A, which results in a first order reaction.

This has appeared when trying to implement the [Robertson problem](https://github.com/DENG-MIT/CRNN/discussions/16) as a test in #817.

Start with _add_reaction() in overreact/core.py and _dydt() in overreact/simulate.py, especially the construction of M from A. Run the reproduction script for "2 B -> B + C" and inspect the existing tests, including the Robertson problem work referenced in #817. Done means the reported reaction order and rate scaling are second order in B without breaking net concentration changes.

Written by the indexing model from the issue text.