Warning for not meeting integration tolerances, while attempting to solve a stiff problem of odes

As far as i can see from other users' questions, this is a quite common warning, when dealing with stiff problems. However, i encountered specific problems despite using appropriate solvers, such as ode15s or ode23s.

Actually, i am dealing with a curve-fitting problem, which requires solving a system of ordinary differential equations.

Here is pretty much what i have done using either fmincon or fminsearch:

%function defining the system of odes to be solved
function dx = FD (t, x, flag, k)
% values from literature
Hthres1 = 0.0003;
Hthres2 = 0.0003;
Hthresmth = 0.080;
Yce1 = 0.006;
Yce2 = 0.006;
Ybut = 0.002;
Yac = 0.001;
Ymth = 0.001;
b = 0.02;
% rate functions
r1 = (k(1)*x(9)*x(1)*(x(8)-Hthres1))/(Yce1*(k(2)+x(1))*(k(3)+(x(8)-Hthres1)))+(k(4)*x(10)*x(1)*(x(8)-Hthres2))/(Yce2*(k(5)+x(1))*(k(6)+(x(8)-Hthres2)));
r2 = (k(1)*x(9)*x(2)*(x(8)-Hthres1))/(Yce1*(k(7)+x(2))*(k(8)+(x(8)-Hthres1)));
r3 = (k(9)*x(9)*x(3))/(k(9)+x(3));
r4 = (k(11)*x(12)*x(6))/(Ybut*(k(12)+x(6)));
r5 = (k(13)*x(13)*x(7))/(Yac*(k(14)+x(7)));
%ODE definitions
dx = zeros (13,1);
dx(1) = -r1;
dx(2) = r1-r2;
dx(3) = r2-r3;
dx(4) = r3;
dx(5) = (0.25*k(15)*x(11)*(x(8)-Hthresmth))/(Ymth*(k(16)+(x(8)-Hthresmth)))+r5;
dx(6) = -r4;
dx(7) = (2*r4)-r5;
dx(8) = (0.5*r4)-r1-r2-(k(15)*x(11)*(x(8)-Hthresmth))/(Ymth*(k(16)+(x(8)-Hthresmth)));
dx(9) = (k(1)*x(9)*x(1)*(x(8)-Hthres1))/((k(2)+x(1))*(k(3)+(x(8)-Hthres1)))+(k(4)*x(10)*x(1)*(x(8)-Hthres2))/((k(5)+x(1))*(k(6)+(x(8)-Hthres2)))-(b*x(9));
dx(10) = (Yce1*r2)-(b*x(10));
dx(11) = (k(15)*x(11)*(x(8)-Hthresmth))/((k(16)+(x(8)-Hthresmth)))-(b*x(11));
dx(12) =  (Ybut*r4) - (b*x(12));
dx(13) = (Yac*r5) - (b*x(13));

end

%function for the system of odes described above
function Err = ED (p, T, TCE, DCE, VC, ETH, MTH)
%solving the system of ODEs
x0 = [TCE(1); DCE(1); VC(1); ETH(1); MTH(1); 300; 85; 0; 9.16; 8.90; 0.29; 2.09; 1.56];
global x
k = p;
options = odeset ('Abstol', 1e-12, 'Reltol', 1e-5);
[t,x] = ode15s('FD',T, x0, options, k);
% Calculating Err
N = length(T);
sum_TCE = 0;
sum_DCE = 0;
sum_VC = 0;
sum_ETH = 0;
sum_MTH = 0;
for i=1:N
  sum_TCE = sum_TCE + (TCE(i) - x(i,1))^2;
  sum_DCE = sum_DCE + (DCE(i)-x(i,2))^2;
  sum_VC = sum_VC + (VC(i)-x(i,3))^2;
  sum_ETH = sum_ETH + (ETH(i) - x(i,4))^2;
  sum_MTH = sum_MTH + (MTH(i) - x(i,5))^2;
end
Err = sqrt (sum_TCE + sum_DCE + sum_VC + sum_ETH + sum_MTH);
end

%script for the curve-fitting process
k = [0.45; 1.5; 0.01; 3; 1.5; 0.007; 10; 0.05; 1.2; 550; 0.16; 20; 0.1; 750; 0.16; 0.8];
lb = [0.3; 0.1; 0.003; 1.0; 0.1; 0.003; 5; 0.01; 0.2; 350; 0.1; 10; 0.02; 150; 0.07; 0.1];
ub = [1.0; 10.0; 0.1; 6; 10; 0.05; 750; 0.1; 10; 1500; 0.25; 50; 0.2; 1200; 0.50; 3];
global x
p0 = k;
options = optimset ('Display', 'iter', 'MaxFunEvals', 9000, 'TolFun', 1e-11, 'TolX', 1e-11, 'maxiter', 500, 'FinDiffType', 'central', 'Algorithm', 'active-set');
[p,resnorm, residual, exitflag, output, lambda, jacobian] = fmincon ('ED', p0, [], [],[], [], lb, ub, [], options, T, TCE, DCE, VC, ETH, MTH);
k = p;
end

The warnings i encounter (_Warning: Failure at t. Unable to meet integration tolerances without reducing the step size below the smallest value allowed_) begin when i change the two boundaries ub and lb, in order to obtain a better fit to my experimental data.

I tried an unconstrained search for the minimum (using fminsearch) and when i changed slightly the initial values (e.g. vector k) i also encountered the same warnings.

Is this problem so "sensitive" to the initial values and the upper and lower boundaries? Or is there something i have written improperly? I'm not considered an experienced Matlab user, so probably there's something beyond my knowledge here and i would really be grateful for any thoughts.

Thank you in advance for your time.