function [sys,x0,str,ts] = sfun_isomerpar(t,x,u,flag,x0,par)



switch flag,
  case 0,
    [sys,x0,str,ts]=mdlInitializeSizes(x0);
  case 1,
    sys=mdlDerivatives(t,x,u);
  case 3,
    sys=mdlOutputs(t,x,u);
  case {2,4,9 }
     sys = []; % Unused flags
  otherwise
    error(['Unhandled flag = ',num2str(flag)]);
end

function [sys,x0,str,ts]=mdlInitializeSizes(x0)
sizes = simsizes;
sizes.NumContStates  = 4;
sizes.NumDiscStates  = 0;
sizes.NumOutputs     = 4;
sizes.NumInputs      = 4;
sizes.DirFeedthrough = 0;
sizes.NumSampleTimes = 1;   % at least one sample time is needed
sys = simsizes(sizes);

%fixed parameters

% initialize the initial conditions
x0=x0; % xinit are the initial values that are sent to the S-function, see first line of the S-function

% str is always an empty matrix
str = [];
% initialize the array of sample times
ts  = [0 0];

function sys=mdlDerivatives(t,x,u)
% Assign values to the parameters, par is a vector with parameter values
% that is sent to the S-function
    umaxl = 0.0887;
    ksl = 0.0194 ;
    kil = 6.502 ;
    umaxz = 0.2016;
    ksz = 4.993;
    kd = 0.0094;
    yl = 0.4106;
    yz= 0.6418;
    alpha=0.0525;
    betha=0.0026;

    par=[umaxl ksl kil umaxz ksz kd yl yz alpha betha];
    umaxl = par(1);
    ksl = par(2);
    kil = par(3);
    umaxz = par(4);
    ksz = par(5);
    kd = par(6);
    yl = par(7);
    yz= par(8);
    alpha= par(9);
    betha= par(10);


uL=umaxl*((x(1)/(ksl+x(1)))*(kil/(kil+x(2))));
uZ=umaxz*(x(2)/(ksz+x(2)));

if t >= 24;
      rp = alpha*uL+betha;
    elseif t <= 24;
       rp = 0;
end



x(4)=rp*x(3);
  
dxdt(1)=-(uL*x(3))/yl;
dxdt(2)=-(uZ*x(3))/yz;
dxdt(3)=uL*x(3)+uZ*x(3)-kd*x(3);
sys = [dxdt(1) dxdt(2) dxdt(3) x(4)];




function sys=mdlOutputs(t,x,u)
% Calculate the outputs of the model block
sys = [x(1) x(2) x(3) x(4)];