Why does defining a system of ODEs a different, but similar way yield very different solutions?

I have a system of ODEs that I am trying solve, however the correct solutions are not being returned (I have another script to which I'm comparing the produced values). As I am trying to figure out why, I defined the function of ODEs as a column vector, instead of how it was originally defined as row vectors, and got a different answer.

My original code produces values closest in magnitude to what I want and is defined as follows:

Enfcn = @(t) Ei + v * t;
kred  = @(t) 225000 * exp(-8) * exp(-(a * f) * (Enfcn(t) - Eo));
kox   = @(t) 225000 * exp(-8) * exp(((1-a) * f) * (Enfcn(t) - Eo));

C0 = [0.1, 0, 0, 0, 0, 1e6, 0];
tspan = 0:0.01:15;

[Time, Concentrations] = ode15s(@(t, C) mechanism(t, C, kred, kox), tspan, C0);

function dC = mechanism(t, C, kred, kox)

%Rate Constants (ish)
k1 = 130;
k2 = 1;
pKa = 2.5;

%Protonation Rates
ka = 1e1;
Ka = 10^(-pKa) * 1e9;
kb = Ka * ka;

    % Variables
        CNiII   = C(1);
        CNiIISH = C(2);
        CNiISH  = C(3);
        CNiIIIH = C(4);
        CNiIIH  = C(5);
        Ce      = C(6); %I know this is unused, but it doesn't hurt having it here
        CH      = C(7);
        CH2     = C(8); %Same with this

     % Rate Laws / ODEs
        dNiII = -ka*CNiII*CH + kb*CNiIISH + k2*CNiIIH*CH;
        dNiIISH = -kb*CNiIISH + ka*CNiII*CH - kred(t)*CNiIISH + kox(t)*CNiISH;
        dNiISH = kred(t)*CNiIISH - kox(t)*CNiISH - k1*CNiISH;
        dNiIIIH = -kred(t)*(CNiIISH + CNiIIIH) + kox(t)*(CNiISH + CNiIIH);
        dNiIIH = k1*CNiISH - kred(t)*CNiIIIH + kox(t)*CNiIIH;
        de = kred(t)*CNiIIIH - kox(t)*CNiIIH - k2*CNiIIH*CH;
        dH = -k2*CNiIIH*CH - ka*CNiII*CH + kb*CNiIISH;
        dH2 = k2*CNiIIH*CH;

     % Assign Output Variables
        dC(1,:) = dNiII;
        dC(2,:) = dNiIISH;
        dC(3,:) = dNiISH;
        dC(4,:) = dNiIIIH;
        dC(5,:) = dNiIIH;
        dC(6,:) = de*(96485.33289/1e9); %This needs to multiplied at some point by this conversion, I do it when I plot, I just added it here
        dC(7,:) = dH;
        dC(8,:) = dH2;

When changing it to this, I gain the plot shape that I want, but not the correct values (in fact this one is several orders of magnitude off), and the integration tolerance of ode15s can't be met even if I set RelTol and AbsTol to values from 1e-12~1e-20:

        dC = [-ka*CNiII*CH + kb*CNiIISH + k2*CNiIIH*CH;
        -kb*CNiIISH + ka*CNiII*CH - kred(t)*CNiIISH + kox(t)*CNiISH;
        kred(t)*CNiIISH - kox(t)*CNiISH - k1*CNiISH;
        -kred(t)*(CNiIISH + CNiIIIH) + kox(t)*(CNiISH + CNiIIH);
        k1*CNiISH - kred(t)*CNiIIIH + kox(t)*CNiIIH;
        kred(t)*CNiIIIH - kox(t)*CNiIIH - k2*CNiIIH*CH;
        -k2*CNiIIH*CH - ka*CNiII*CH + kb*CNiIISH;
        k2*CNiIIH*CH];

I am confused as to why these give different results when from my understanding they were synonymous, but I am also still new to this. Any insight would be very appreciated. Thanks in advance.

Answer 1

[James Tursa]

Do you simply need to apply that factor in the 6th spot?

(kred(t)*CNiIIIH - kox(t)*CNiIIH - k2*CNiIIH*CH)*(96485.33289/1e9);

[Riley Stein]

I figured out the problem because of your comment. The original function can be in any order because I later define what is exaclty what, whereas the column vector is sensitive to exactly where the ODE is placed, which makes sense. I'm sorry for such a senseless error. I was too excited about getting the plot shape correct that I didn't even consider the actual variable naming. Although I am still perplexed by why the ODEs are not being solved the same.