Fix calculating jacobian from mix of rate rules and reactions.

[luciansmith]

The following two scripts work fine:

r = te.loada ('''
      species x
      x' = -k1*x
      k1 = 0.1; x = 1
''')

print(r.getFullJacobian())

gives you:

        x
x [[ -0.1]]

and:

r = te.loada ('''
      S1 -> S2; k3*S1
      S1 = 10
      k3=0.5
''')

print(r.getFullJacobian())

gives you:

        S1, S2
S1 [[ -0.5,  0],
S2  [  0.5,  0]]

but if you try to combine them:

r = te.loada ('''
      species x
      x' = -k1*x
      k1 = 0.1; x = 1
      S1 -> S2; k3*S1
      S1 = 10
      k3=0.5
''')

print(r.getFullJacobian())

you get:

        S1, S2, x
S1 [[    0,  0, 0],
S2  [ -0.5,  0, 0],
x   [  0.5,  0, 0]]

which gets the row values wrong for the 0.5, -0.5, and misses out [x,x] entirely.

If you make 'x' a boundary species, it falls back out of the Jacobian and you get the correct values again, but only for S1/S2.