How can I improve my non linear multi dimensional regression? 8 parameters and 8 sets of data, 4 variables

I have this system of equations that I need to fit 8 parameters with 8/maybe 9 sets of experimental data:

My response variable is [DCF], I have three independent variables: [O_3], [FeOOH] and time. The variable [OH] is a function that depends of [DCF] and must be evaluated as third degree polynomial and then inserted into the differential equation. Since I didn't know how to create a code that selectively chose the correct root of the third degree polynomial during the regression I instead derived the equation implicitly and got this system of differential equations:

I have 8 parameters that I have represented in the following way:

I followed the advice provided by user @Star Strider and managed to fit my experimental data, 8 sets of data to my set of differential equations. However, every set of data has differente values for [FeOOH] and [O_3] and so, they cannot be bundled together but rather individually. Attatched is my code. When I began with the first set of data, it fitted perfectly:

But now that I've included my 8 sets of data I cannot achieve that same curve for this data:

Is there any way to improve the efficiency or accuracy of lsqcurvefit to better fit my data? I've already tried ramping up the tolerances using optimoptions and odeset for my ODE's and for the parameter estimation. I also tried setting the initial condition of dx/dt as parameters to be fitted (because dx/dt is a theoretical value anyways) and altough it fits nicely in every curve, it's not the case for Curve 1.

Thanks in advance for any help