Find the steady states of a system.

I am building a SimBiology Model. I want to simulate from randomly generated initial conditions numerous times and find the steady state of a particular specie (PLS). MY current code does this in a graphical sense but I need the actual numbers in a vector (for example) for subsequent use. Does anybody know how to do this? My current code is:

%Random Number Generations and Assignment
No_of_simulations = 10;      %i.e. points
number_of_parameters = 4;   %i.e. dimenions
lb=[1e-3];
ub=[0.1];
x = lhsdesign(number_of_parameters ,No_of_simulations);
D = bsxfun(@plus,lb,bsxfun(@times,x,(ub-lb)));

for i=x,
    n1 = i(1);
    n2 = i(2);
    n3 = i(3);
    n4 = i(4);   

%Model generation. 

      Mobj = sbiomodel('u-PA_parameter_generation');
      comp_obj = addcompartment(Mobj, 'plasma');
      %Reaction 1
      Robj1 = addreaction(Mobj, 'Pro-u-PA + PLG -> PLS + Pro-u-PA');
      Kobj1= addkineticlaw(Robj1, 'MassAction');
      Pobj1 = addparameter(Kobj1, 'keff_zymogen',0.035); 
      set(Kobj1, 'ParameterVariableNames','keff_zymogen');

      %Reaction 2
      Robj2 = addreaction(Mobj, 'PLS + Pro-u-PA -> PLS + u-PA');
      Kobj2= addkineticlaw(Robj2, 'MassAction');
      Pobj2 = addparameter(Kobj2, 'keff_PLS',40);
      set(Kobj2, 'ParameterVariableNames','keff_PLS');

      %Reaction 3
      Robj3 = addreaction(Mobj, 'u-PA + PLG -> u-PA + PLS');
      Kobj3= addkineticlaw(Robj3, 'MassAction');
      Pobj3 = addparameter(Kobj3, 'keff_pos',0.9);
      set(Kobj3, 'ParameterVariableNames','keff_pos');

      %Reaction 4
      Robj4 = addreaction(Mobj, 'Pro-u-PA -> null');
      Kobj4= addkineticlaw(Robj4, 'MassAction');
      Pobj4 = addparameter(Kobj4, 'u1',0.084);
      set(Kobj4, 'ParameterVariableNames','u1');

      %Reaction 5
      Robj5 = addreaction(Mobj, 'PLG -> null');
      Kobj5= addkineticlaw(Robj5, 'MassAction');
      Pobj5 = addparameter(Kobj5, 'u2',0.032);
      set(Kobj5, 'ParameterVariableNames','u2');

      %Reaction 6
      Robj6 = addreaction(Mobj, 'PLS -> null');
      Kobj6= addkineticlaw(Robj6, 'MassAction');
      Pobj6 = addparameter(Kobj6, 'u1',0.084); %Same as R4
      set(Kobj6, 'ParameterVariableNames','u1');

      %Reaction 7
      Robj7 = addreaction(Mobj, 'u-PA -> null');
      Kobj7= addkineticlaw(Robj7, 'MassAction');
      Pobj7 = addparameter(Kobj7, 'u1',0.084); %Same as R4
      set(Kobj7, 'ParameterVariableNames','u1');

      %Reaction 8
      Robj8 = addreaction(Mobj, 'null -> Pro-u-PA');
      Kobj8= addkineticlaw(Robj8, 'MassAction');
      Pobj8 = addparameter(Kobj8, 'a1',0.0032);
      set(Kobj8, 'ParameterVariableNames','a1');

      %Reaction 9
      Robj9 = addreaction(Mobj, 'null -> PLG');
      Kobj9= addkineticlaw(Robj9, 'MassAction');
      Pobj9 = addparameter(Kobj9, 'a2',0.01);
      set(Kobj9, 'ParameterVariableNames','a2');

%setting species concentrations
    Sobj1 = sbioselect(Mobj,'Type','species','Name','Pro-u-PA');
    set(Sobj1, 'InitialAmount',n1)

      Sobj2 = sbioselect(Mobj,'Type','species','Name','PLG');
      set(Sobj2, 'InitialAmount',n2)

      Sobj3 = sbioselect(Mobj,'Type','species','Name','PLS');
      set(Sobj3, 'InitialAmount',n3)

      Sobj4 = sbioselect(Mobj,'Type','species','Name','u-PA');
      set(Sobj4, 'InitialAmount',n4)

%simulate
    config = getconfigset(Mobj);
    set(config,'StopTime',1000)
    [t_ode, x_ode, names] = sbiosimulate(Mobj); 
    hold on
    figure;
    set(gcf, 'color','white');
    plot(t_ode, x_ode(:,1:end));
    legend(names)
end

Thanks in Advance

Answer 1

[Ciaran]

Basically the answer was to initialize the figure before the for loop

Answer 2

[Ingrid Tigges]

Given that the steady state is defined as dx/dt=0 which is approximately delta x/delta t you can check whether the difference in x of two consecutive time points is 0

diff_x= diff(x_ode);

Due to numeric errors you want to have those differences that are below a certain threshold:

 threshold = 1e-15;
 x_equals_0 = diff_x<threshold;

The indices of those values in x_equals_0 that are 1 belong to points where dx/dt=0.