Stability analysis of a non-linear ODE system

I solved the following ODE system using the code:

syms Sci C Sr Sh R Cf Cp Ce E HR H Sp P k1 k2 k3 k4 k5 k6 k7 k8 k9 k10 k11 k12 k13 k14 k15 k16 p1 p2 p3 mu eta theta alpha CL
F=zeros(13,1);
e1=F(1) == (mu - eta*Sci*C);
e2=F(1) == (theta*Ce - eta*Sci*C);
e3=F(3) == (k1*Sci - k2*Sr);
e4=F(4) == (k1*Sci - k10*Sh);
e5=F(5) == (p1*Sr - k11*R - k14*P*R);
e6=F(6) == (k3*CL*R - k4*Cf);
e7=F(7) == (k4*Cf - k5*Cp + k6*Ce);
e8=F(8) == (k5*Cp - k6*Ce + k9*HR*H - k12*Ce - k7*Ce + k8*E);
e9=F(9) == (k7*Ce - k8*E);
e10=F(10) == (p2*Sh - k15*HR*Ce);
e11=F(11) == (alpha - k9*HR*H);
e12=F(12) == (k1*Sci - k13*Sp);
e13=F(13) == (p3*Sp - k16*P);

[Sci,C,Sr,Sh,R,Cf,Cp,Ce,E,HR,H,Sp,P]=solve(e1,e2,e3,e4,e5,e6,e7,e8,e9,e10,e11,e12,e13,Sci,C,Sr,Sh,R,Cf,Cp,Ce,E,HR,H,Sp,P)

Sci = 

C = 

Sr = 

Sh = 

R = 

Cf = 

Cp = 

Ce = 

E = 

HR = 

H = 

Sp = 

P = 

with stability points for C_e = mu/theta, C_p = (k6*mu - alpha*theta + k12*mu)/(k5*theta), C_f = -(alpha*theta - k12*mu)/(k4*theta) and E = (k7*mu)/(k8*theta). All other variables are dependent on an input variable labelled CL. I want to do a stability analysis around these points.

I've determined the Jacobian using the code:

syms Sci C Sr Sh R Cf Cp Ce E HR H Sp P k1 k2 k3 k4 k5 k6 k7 k8 k9 k10 k11 k12 k13 k14 k15 k16 p1 p2 p3 mu eta theta alpha CL

f = [mu - eta*Sci*C,theta*Ce - eta*Sci*C,k1*Sci - k2*Sr,k1*Sci - k10*Sh,p1*Sr - k11*R - k14*P*R,k3*CL*R - k4*Cf,k4*Cf - k5*Cp + k6*Ce,k5*Cp - k6*Ce + k9*HR*H - k12*Ce - k7*Ce + k8*E,k7*Ce - k8*E,p2*Sh - k15*HR*Ce,alpha - k9*HR*H,k1*Sci - k13*Sp,p3*Sp - k16*P];

J = jacobian(f, [Sci,C,Sr,Sh,R,Cf,Cp,Ce,E,HR,H,Sp,P]);

J = simplify(J)

J = 

I now have to determine the eigenvalues by first substituting equilibrium point values for the variables in the Jacobian with the steady state values ( in terms of CL) generated in the code above. I'm not sure how to do this in Matlab?

Answer 1

[Torsten]

You won't succeed in this generality. To determine the eigenvalues, MATLAB had to solve for the roots of a polynomial of degree 13 with symbolic coefficients. This is in general only possible for polynomials up to degree 4. So you have to give values to the parameters of your function, I guess.

syms Sci C Sr Sh R Cf Cp Ce E HR H Sp P k1 k2 k3 k4 k5 k6 k7 k8 k9 k10 k11 k12 k13 k14 k15 k16 p1 p2 p3 mu eta theta alpha CL
F=zeros(13,1);
e1=F(1) == (mu - eta*Sci*C);
e2=F(1) == (theta*Ce - eta*Sci*C);
e3=F(3) == (k1*Sci - k2*Sr);
e4=F(4) == (k1*Sci - k10*Sh);
e5=F(5) == (p1*Sr - k11*R - k14*P*R);
e6=F(6) == (k3*CL*R - k4*Cf);
e7=F(7) == (k4*Cf - k5*Cp + k6*Ce);
e8=F(8) == (k5*Cp - k6*Ce + k9*HR*H - k12*Ce - k7*Ce + k8*E);
e9=F(9) == (k7*Ce - k8*E);
e10=F(10) == (p2*Sh - k15*HR*Ce);
e11=F(11) == (alpha - k9*HR*H);
e12=F(12) == (k1*Sci - k13*Sp);
e13=F(13) == (p3*Sp - k16*P);

[Scieq,Ceq,Sreq,Sheq,Req,Cfeq,Cpeq,Ceeq,Eeq,HReq,Heq,Speq,Peq]=solve([e1,e2,e3,e4,e5,e6,e7,e8,e9,e10,e11,e12,e13],[Sci,C,Sr,Sh,R,Cf,Cp,Ce,E,HR,H,Sp,P]);

f = [rhs(e1),rhs(e2),rhs(e3),rhs(e4),rhs(e5),rhs(e6),rhs(e7),rhs(e8),rhs(e9),rhs(e10),rhs(e11),rhs(e12),rhs(e13)];
%f = [mu - eta*Sci*C,theta*Ce - eta*Sci*C,k1*Sci - k2*Sr,k1*Sci - k10*Sh,p1*Sr - k11*R - k14*P*R,k3*CL*R - k4*Cf,k4*Cf - k5*Cp + k6*Ce,k5*Cp - k6*Ce + k9*HR*H - k12*Ce - k7*Ce + k8*E,k7*Ce - k8*E,p2*Sh - k15*HR*Ce,alpha - k9*HR*H,k1*Sci - k13*Sp,p3*Sp - k16*P];

J = jacobian(f, [Sci,C,Sr,Sh,R,Cf,Cp,Ce,E,HR,H,Sp,P]);

J = simplify(J)

Jeq = subs(J,[Sci,C,Sr,Sh,R,Cf,Cp,Ce,E,HR,H,Sp,P],[Scieq,Ceq,Sreq,Sheq,Req,Cfeq,Cpeq,Ceeq,Eeq,HReq,Heq,Speq,Peq])
eig(Jeq)

[Ron_S]

Thank you!

[Ron_S]

Further to my initial question, I have now given specific values for my parameters and managed to find the eigen values by finding the roots of the characteristic polynomial generated.

I would now like to run the code for a collection of randomly chosen parameter values (say a hundred or so different parameter sets) within a for loop, and check that the eigenvalues have negative real part in each case. I just don't know how to go about doing that or if it is possible?

Thanks in advance,

Ron

%Using the parameter values, with CL = 0.5
k1 = 3;
k2 = 2;
k3 = 5;
k4 = 4;
k5 = 5;
k6 = 1;
k7 = 4;
k8 = 3;
k9 = 1;
k10 = 1.2;
k11 = 1;
k12 = 9;
k13 = 1;
k14 = 1;
k15 = 1;
k16 = 1;
p1 = 30;
p2 = 4;
p3 = 1.5;
mu = 25;
eta = 10;
alpha = 10;
theta = 10;
CL = 0.5;

syms Sci C Sr Sh R Cf Cp Ce E HR H Sp P

F=zeros(13,1);
e1=F(1) == (mu - eta*Sci*C);
e2=F(1) == (theta*Ce - eta*Sci*C);
e3=F(3) == (k1*Sci - k2*Sr);
e4=F(4) == (k1*Sci - k10*Sh);
e5=F(5) == (p1*Sr - k11*R - k14*P*R);
e6=F(6) == (k3*CL*R - k4*Cf);
e7=F(7) == (k4*Cf - k5*Cp + k6*Ce);
e8=F(8) == (k5*Cp - k6*Ce + k9*HR*H - k12*Ce - k7*Ce + k8*E);
e9=F(9) == (k7*Ce - k8*E);
e10=F(10) == (p2*Sh - k15*HR*Ce);
e11=F(11) == (alpha - k9*HR*H);
e12=F(12) == (k1*Sci - k13*Sp);
e13=F(13) == (p3*Sp - k16*P);

%Now solve to find the steady state values in terms of (A, B, C, D, Cl and the parameters as listed), then determine the Jacobian J:

[Scieq,Ceq,Sreq,Sheq,Req,Cfeq,Cpeq,Ceeq,Eeq,HReq,Heq,Speq,Peq]=solve([e1,e2,e3,e4,e5,e6,e7,e8,e9,e10,e11,e12,e13],[Sci,C,Sr,Sh,R,Cf,Cp,Ce,E,HR,H,Sp,P]);
f = [rhs(e1),rhs(e2),rhs(e3),rhs(e4),rhs(e5),rhs(e6),rhs(e7),rhs(e8),rhs(e9),rhs(e10),rhs(e11),rhs(e12),rhs(e13)];
J = jacobian(f, [Sci,C,Sr,Sh,R,Cf,Cp,Ce,E,HR,H,Sp,P]);
J = simplify(J)

J = 

%Substitute the steady state values obtained into the Jacobian:

Jeq = subs(J,[Sci,C,Sr,Sh,R,Cf,Cp,Ce,E,HR,H,Sp,P],[Scieq,Ceq,Sreq,Sheq,Req,Cfeq,Cpeq,Ceeq,Eeq,HReq,Heq,Speq,Peq])

Jeq = 

%Find the eigenvalue matrix:
eig(Jeq)

ans = 

%Matlab does not solve the eigenvalues of symbolic matrices greater that 4 x 4. Instead, it gives it as the characteristic polynomial, in this case, a 13 degree polynomial, (z^13...).
%We can find the roots of this polynomial as follow:
%List the coefficients from the characteristic polynomial sigma_1 as a vector and then find the roots:
p = [1 (6809/45) 7658861/1620 17552999/270 270831481/540 3922302799/1620 1253029309/162 27428625103/1620 21201142273/810 2376322271/81 1924251074/81 347643560/27 11572400/3 400000]

p = 1×14

1.0e+07 * 0.0000 0.0000 0.0005 0.0065 0.0502 0.2421 0.7735 1.6931 2.6174 2.9337 2.3756 1.2876 0.3857 0.0400

r = roots(p)

r = 
   1.0e+02 *

  -1.1472 + 0.0000i
  -0.1465 + 0.0000i
  -0.0550 + 0.0000i
  -0.0337 + 0.0182i
  -0.0337 - 0.0182i
  -0.0307 + 0.0000i
  -0.0244 + 0.0000i
  -0.0022 + 0.0113i
  -0.0022 - 0.0113i
  -0.0152 + 0.0000i
  -0.0113 + 0.0000i
  -0.0091 + 0.0000i
  -0.0019 + 0.0000i

Answer 2

[Torsten]

syms Sci C Sr Sh R Cf Cp Ce E HR H Sp P k1 k2 k3 k4 k5 k6 k7 k8 k9 k10 k11 k12 k13 k14 k15 k16 p1 p2 p3 mu eta alpha theta CL
F=zeros(13,1);
e1=F(1) == (mu - eta*Sci*C);
e2=F(1) == (theta*Ce - eta*Sci*C);
e3=F(3) == (k1*Sci - k2*Sr);
e4=F(4) == (k1*Sci - k10*Sh);
e5=F(5) == (p1*Sr - k11*R - k14*P*R);
e6=F(6) == (k3*CL*R - k4*Cf);
e7=F(7) == (k4*Cf - k5*Cp + k6*Ce);
e8=F(8) == (k5*Cp - k6*Ce + k9*HR*H - k12*Ce - k7*Ce + k8*E);
e9=F(9) == (k7*Ce - k8*E);
e10=F(10) == (p2*Sh - k15*HR*Ce);
e11=F(11) == (alpha - k9*HR*H);
e12=F(12) == (k1*Sci - k13*Sp);
e13=F(13) == (p3*Sp - k16*P);
%Now solve to find the steady state values in terms of (A, B, C, D, Cl and the parameters as listed), then determine the Jacobian J:
[Scieq,Ceq,Sreq,Sheq,Req,Cfeq,Cpeq,Ceeq,Eeq,HReq,Heq,Speq,Peq]=solve([e1,e2,e3,e4,e5,e6,e7,e8,e9,e10,e11,e12,e13],[Sci,C,Sr,Sh,R,Cf,Cp,Ce,E,HR,H,Sp,P]);
f = [rhs(e1),rhs(e2),rhs(e3),rhs(e4),rhs(e5),rhs(e6),rhs(e7),rhs(e8),rhs(e9),rhs(e10),rhs(e11),rhs(e12),rhs(e13)];
J = jacobian(f, [Sci,C,Sr,Sh,R,Cf,Cp,Ce,E,HR,H,Sp,P]);
J = simplify(J) 

J = 

%Substitute the steady state values obtained into the Jacobian:
Jeq = subs(J,[Sci,C,Sr,Sh,R,Cf,Cp,Ce,E,HR,H,Sp,P],[Scieq,Ceq,Sreq,Sheq,Req,Cfeq,Cpeq,Ceeq,Eeq,HReq,Heq,Speq,Peq])

Jeq = 

%Lower (lb) and upper (ub) bounds for the free parameters
lbk1 = 3-0.5*3;
lbk2 = 2-0.5*2;
lbk3 = 5-0.5*5;
lbk4 = 4-0.5*4;
lbk5 = 5-0.5*5;
lbk6 = 1-0.5*1;
lbk7 = 4-0.5*4;
lbk8 = 3-0.5*3;
lbk9 = 1-0.5*1;
lbk10 = 1.2-0.5*1.2;
lbk11 = 1-0.5*1;
lbk12 = 9-0.5*9;
lbk13 = 1-0.5*1;
lbk14 = 1-0.5*1;
lbk15 = 1-0.5*1;
lbk16 = 1-0.5*1;
lbp1 = 30-0.5*30;
lbp2 = 4-0.5*4;
lbp3 = 1.5-0.5*1.5;
lbmu = 25-0.5*25;
lbeta = 10-0.5*10;
lbalpha = 10-0.5*10;
lbtheta = 10-0.5*10;
lbCL = 0.5-0.5*0.5;

ubk1 = 3+0.5*3;
ubk2 = 2+0.5*2;
ubk3 = 5+0.5*5;
ubk4 = 4+0.5*4;
ubk5 = 5+0.5*5;
ubk6 = 1+0.5*1;
ubk7 = 4+0.5*4;
ubk8 = 3+0.5*3;
ubk9 = 1+0.5*1;
ubk10 = 1.2+0.5*1.2;
ubk11 = 1+0.5*1;
ubk12 = 9+0.5*9;
ubk13 = 1+0.5*1;
ubk14 = 1+0.5*1;
ubk15 = 1+0.5*1;
ubk16 = 1+0.5*1;
ubp1 = 30+0.5*30;
ubp2 = 4+0.5*4;
ubp3 = 1.5+0.5*1.5;
ubmu = 25+0.5*25;
ubeta = 10+0.5*10;
ubalpha = 10+0.5*10;
ubtheta = 10+0.5*10;
ubCL = 0.5+0.5*0.5;

%Number of trials
n = 100;
found = zeros(n,1);

% Make a Monte Carlo simulation
for i = 1:n
  % Random values for the parameters uniformly distributed between their
  % lower and upper bounds
  k1num(i) = lbk1 + rand()*(ubk1-lbk1); 
  k2num(i) = lbk2 + rand()*(ubk2-lbk2);
  k3num(i) = lbk3 + rand()*(ubk3-lbk3); 
  k4num(i) = lbk4 + rand()*(ubk4-lbk4);
  k5num(i) = lbk5 + rand()*(ubk5-lbk5); 
  k6num(i) = lbk6 + rand()*(ubk6-lbk6);
  k7num(i) = lbk7 + rand()*(ubk7-lbk7); 
  k8num(i) = lbk8 + rand()*(ubk8-lbk8);
  k9num(i) = lbk9 + rand()*(ubk9-lbk9); 
  k10num(i) = lbk10 + rand()*(ubk10-lbk10);
  k11num(i) = lbk11 + rand()*(ubk11-lbk11); 
  k12num(i) = lbk12 + rand()*(ubk12-lbk12);
  k13num(i) = lbk13 + rand()*(ubk13-lbk13); 
  k14num(i) = lbk14 + rand()*(ubk14-lbk14);
  k15num(i) = lbk15 + rand()*(ubk15-lbk15); 
  k16num(i) = lbk16 + rand()*(ubk16-lbk16);
  p1num(i) = lbp1 + rand()*(ubp1-lbp1); 
  p2num(i) = lbp2 + rand()*(ubp2-lbp2);
  p3num(i) = lbp3 + rand()*(ubp3-lbp3);
  munum(i) = lbmu + rand()*(ubmu-lbmu); 
  etanum(i) = lbeta + rand()*(ubeta-lbeta);
  alphanum(i) = lbalpha + rand()*(ubalpha-lbalpha); 
  thetanum(i) = lbtheta + rand()*(ubtheta-lbtheta);
  CLnum(i) = lbCL + rand()*(ubCL-lbCL);
  Jeqnum = subs(Jeq,[k1 k2 k3 k4 k5 k6 k7 k8 k9 k10 k11 k12 k13 k14 k15 k16 p1 p2 p3 mu eta alpha theta CL],[k1num(i) k2num(i) k3num(i) k4num(i) k5num(i) k6num(i) k7num(i) k8num(i) k9num(i) k10num(i) k11num(i) k12num(i) k13num(i) k14num(i) k15num(i) k16num(i) p1num(i) p2num(i) p3num(i) munum(i) etanum(i) alphanum(i) thetanum(i) CLnum(i)]);
  % Compute eigenvalues
  v{i} = double(eig(Jeqnum));
  % Check if all real parts are <= 0
  if any(real(v{i})>0) 
    found(i) = 1;
  end
end
% List the number of trials where eigenvalues with real part > 0 occured
sum(found)

ans = 41

[Ron_S]

Thank you so much Torsten, this makes sense!

[Ron_S]

I have another question to ask please :-)

Would it be possible to run a trial and write some code to extract a specific set of parameter values for which the eigenvalues have real parts > 0? I have no idea how to write the code for this, if possible.

Thank you

[Torsten]

I don't understand exactly what you mean. Could you explain ?

[Sam Chak]

@Ron_S: Would it be possible to run a trial and write some code to extract a specific (what?) set of parameter values (from what object or structure?) for which the eigenvalues have real parts > 0?

Are you referring the Jacobian matrix of the set of nonlinear ODEs?

If so, what do the eigenvalues of the Jacobian matrix with positive real parts imply?

And, how about for the rest of the eigenvalues with zero and negative real parts?

[Ron_S]

@Torsten, @Sam Chak, my apologies, I'll restate my question.

When running the current code with random parameter values, the output lists the number of trials where eigenvalues with real part > 0 occurred. For example, out of 100 trials, there could be 2 trials where eigenvalues with real part > 0 occurred.

Is it possible to find the values of the parameters that was substituted into the Jacobian in those two trials specifically?

The eigenvalues of the Jacobian matrix with positive real parts would imply an unstable steady state of the system. I would like to use those values in another simulation.

[Torsten]

found(i) and parameters(i) correspond. So if found(i) == 1, just output k1num(i), k2num(i),... to get a set of parameters for which at least one eigenvalue of the Jacobian is > 0.

[Ron_S]

Dear Torsten,

I'm not 100% sure the code that I now have generates a set of parameters for which at least one eigenvalue of the Jacobian is > 0? I've added the code in lthe last line.

When I check the specific k1num(i) value in the 1x100 double in the workspace that corresponds with the trial number that has a eigenvalue >0 it does not correspond? (Screenshot attached)

Thanks for all your help!

clear
syms Sci C Sr Sh R Cf Cp Ce E HR H Sp P k1 k2 k3 k4 k5 k6 k7 k8 k9 k10 k11 k12 k13 k14 k15 k16 p1 p2 p3 mu eta alpha theta CL
F=zeros(13,1);
e1=F(1) == (mu - eta*Sci*C);
e2=F(1) == (theta*Ce - eta*Sci*C);
e3=F(3) == (k1*Sci - k2*Sr);
e4=F(4) == (k1*Sci - k10*Sh);
e5=F(5) == (p1*Sr - k11*R - k14*P*R);
e6=F(6) == (k3*CL*R - k4*Cf);
e7=F(7) == (k4*Cf - k5*Cp + k6*Ce);
e8=F(8) == (k5*Cp - k6*Ce + k9*HR*H - k12*Ce - k7*Ce + k8*E);
e9=F(9) == (k7*Ce - k8*E);
e10=F(10) == (p2*Sh - k15*HR*Ce);
e11=F(11) == (alpha - k9*HR*H);
e12=F(12) == (k1*Sci - k13*Sp);
e13=F(13) == (p3*Sp - k16*P);
%Now solve to find the steady state values in terms of (A, B, C, D, Cl and the parameters as listed), then determine the Jacobian J:
[Scieq,Ceq,Sreq,Sheq,Req,Cfeq,Cpeq,Ceeq,Eeq,HReq,Heq,Speq,Peq]=solve([e1,e2,e3,e4,e5,e6,e7,e8,e9,e10,e11,e12,e13],[Sci,C,Sr,Sh,R,Cf,Cp,Ce,E,HR,H,Sp,P]);
f = [rhs(e1),rhs(e2),rhs(e3),rhs(e4),rhs(e5),rhs(e6),rhs(e7),rhs(e8),rhs(e9),rhs(e10),rhs(e11),rhs(e12),rhs(e13)];
J = jacobian(f, [Sci,C,Sr,Sh,R,Cf,Cp,Ce,E,HR,H,Sp,P]);
J = simplify(J) 
%Substitute the steady state values obtained into the Jacobian:
Jeq = subs(J,[Sci,C,Sr,Sh,R,Cf,Cp,Ce,E,HR,H,Sp,P],[Scieq,Ceq,Sreq,Sheq,Req,Cfeq,Cpeq,Ceeq,Eeq,HReq,Heq,Speq,Peq])
%Lower (lb) and upper (ub) bounds for the free parameters
lbk1 = 3-0.5*3;
lbk2 = 2-0.5*2;
lbk3 = 5-0.5*5;
lbk4 = 4-0.5*4;
lbk5 = 5-0.5*5;
lbk6 = 1-0.5*1;
lbk7 = 4-0.5*4;
lbk8 = 3-0.5*3;
lbk9 = 1-0.5*1;
lbk10 = 1.2-0.5*1.2;
lbk11 = 1-0.5*1;
lbk12 = 9-0.5*9;
lbk13 = 1-0.5*1;
lbk14 = 1-0.5*1;
lbk15 = 1-0.5*1;
lbk16 = 1-0.5*1;
lbp1 = 30-0.5*30;
lbp2 = 4-0.5*4;
lbp3 = 1.5-0.5*1.5;
lbmu = 25-0.5*25;
lbeta = 10-0.5*10;
lbalpha = 10-0.5*10;
lbtheta = 6-0.5*6;
lbCL = 9-0.5*9;

ubk1 = 3+0.5*3;
ubk2 = 2+0.5*2;
ubk3 = 5+0.5*5;
ubk4 = 4+0.5*4;
ubk5 = 5+0.5*5;
ubk6 = 1+0.5*1;
ubk7 = 4+0.5*4;
ubk8 = 3+0.5*3;
ubk9 = 1+0.5*1;
ubk10 = 1.2+0.5*1.2;
ubk11 = 1+0.5*1;
ubk12 = 9+0.5*9;
ubk13 = 1+0.5*1;
ubk14 = 1+0.5*1;
ubk15 = 1+0.5*1;
ubk16 = 1+0.5*1;
ubp1 = 30+0.5*30;
ubp2 = 4+0.5*4;
ubp3 = 1.5+0.5*1.5;
ubmu = 25+0.5*25;
ubeta = 10+0.5*10;
ubalpha = 10+0.5*10;
ubtheta = 6+0.5*6;
ubCL = 9+0.5*9;

%Number of trials
n = 100;
found = zeros(n,1);

% Make a Monte Carlo simulation
for i = 1:n
  % Random values for the parameters uniformly distributed between their
  % lower and upper bounds
  k1num(i) = lbk1 + rand()*(ubk1-lbk1); 
  k2num(i) = lbk2 + rand()*(ubk2-lbk2);
  k3num(i) = lbk3 + rand()*(ubk3-lbk3); 
  k4num(i) = lbk4 + rand()*(ubk4-lbk4);
  k5num(i) = lbk5 + rand()*(ubk5-lbk5); 
  k6num(i) = lbk6 + rand()*(ubk6-lbk6);
  k7num(i) = lbk7 + rand()*(ubk7-lbk7); 
  k8num(i) = lbk8 + rand()*(ubk8-lbk8);
  k9num(i) = lbk9 + rand()*(ubk9-lbk9); 
  k10num(i) = lbk10 + rand()*(ubk10-lbk10);
  k11num(i) = lbk11 + rand()*(ubk11-lbk11); 
  k12num(i) = lbk12 + rand()*(ubk12-lbk12);
  k13num(i) = lbk13 + rand()*(ubk13-lbk13); 
  k14num(i) = lbk14 + rand()*(ubk14-lbk14);
  k15num(i) = lbk15 + rand()*(ubk15-lbk15); 
  k16num(i) = lbk16 + rand()*(ubk16-lbk16);
  p1num(i) = lbp1 + rand()*(ubp1-lbp1); 
  p2num(i) = lbp2 + rand()*(ubp2-lbp2);
  p3num(i) = lbp3 + rand()*(ubp3-lbp3);
  munum(i) = lbmu + rand()*(ubmu-lbmu); 
  etanum(i) = lbeta + rand()*(ubeta-lbeta);
  alphanum(i) = lbalpha + rand()*(ubalpha-lbalpha); 
  thetanum(i) = lbtheta + rand()*(ubtheta-lbtheta);
  CLnum(i) = lbCL + rand()*(ubCL-lbCL);
  Jeqnum = subs(Jeq,[k1 k2 k3 k4 k5 k6 k7 k8 k9 k10 k11 k12 k13 k14 k15 k16 p1 p2 p3 mu eta alpha theta CL], ...
      [k1num(i) k2num(i) k3num(i) k4num(i) k5num(i) k6num(i) k7num(i) k8num(i) k9num(i) k10num(i) k11num(i) k12num(i) ...
      k13num(i) k14num(i) k15num(i) k16num(i) p1num(i) p2num(i) p3num(i) munum(i) etanum(i) alphanum(i) thetanum(i) CLnum(i)]);
  % Compute eigenvalues
  v{i} = double(eig(Jeqnum));
  % Check if all real parts are <= 0
  if any(real(v{i})>0) 
    found(i) = 1;
  end
end
% List the number of trials where eigenvalues with real part > 0 occured
sum(found)
% To get a set of parameters for which at least one eigenvalue of the Jacobian is > 0.
k1num(i), k2num(i), k3num(i), k4num(i), k5num(i), k6num(i), k7num(i), k8num(i), k9num(i), k10num(i), k11num(i), k12num(i), k13num(i), k14num(i), k15num(i), k16num(i), p1num(i), p2num(i), p3num(i), munum(i), etanum(i), alphanum(i), thetanum(i), CLnum(i),(found)

[Torsten]

Your output of parameters is after the loop over i, thus for i = n.

If found(n) = 1 by chance, you got a set of parameters for which at least one eigenvalue of the Jacobian is > 0.

But this will usually not be the case.

Try to understand the code and why you have to output the parameters here:

if any(real(v{i})>0) 
    found(i) = 1;
    k1num(i),k2num(i),....
end