function [sys,x0,str,ts] = rea7m1(t,x,u,flag)
% Non-linear deviation model of the CSTR;
% 3 equations, Reaction of the type A ----> B;

%   k0 = a1       [1/min];
%   g = a2        [K];
%   ro = a3       [kg/m3];
%   cp = a4       [kJ/kg.K];
%   V = a5        [m3];
%   Hr = a6       [kJ/kmol];
%   F = a7        [m2];
%   al = a8       [kJ/m2.min.K];
%   roc = a9      [kg/m3];
%   cpc = a10     [kJ/kg.K];
%   Vc = a11      [m3];
%   qs = a12      [m3/min];
%   Thvs = a13    [K];
%   cavs = a14    [kmol/m3];
%   Thcvs = a15   [K];
%   qcs = a16     [m3/min];
%   cas = a17     [kmol/m3];
%   Ths = a18     [K];
%   Thcs = a19    [K];

%   x(1) - concentration of A, dimensionless, deviation;
%   x(2) - temperature of the mixture, dimensionless, deviation; 
%   x(3) - jacket temperature, dimensionless, deviation;
%   u(1) - cooling flowrate, dimensionless, deviation;
%   u=qc/qcs ;
%   u(1)=<-0.2,   0.2> = cca<-0.001m3/min,0.001m3/min>;

  switch flag,
   case 0     
    [sys,x0,str,ts] = mdlInitializeSizes;
   case 1,
    sys = mdlDerivatives(t,x,u);
   case 3,
    sys = mdlOutputs(t,x,u);
   case 9                                                
    sys = [];
  end
  
function [sys,x0,str,ts] = mdlInitializeSizes()
  sizes = simsizes;
  
  sizes.NumContStates  = 3;
  sizes.NumDiscStates  = 0;
  sizes.NumOutputs     = 6;
  sizes.NumInputs      = 1;
  sizes.DirFeedthrough = 0;
  sizes.NumSampleTimes = 1;
  
  sys = simsizes(sizes);
  x0  = [0; 0; 0];
  str = [];
  ts  = [0 0];
  
function sys = mdlDerivatives(t,x,u)
% Calculation of f1c, f2c,...,f7c in the steady-state;
  a1=10000000000000;
  a2=11078;
  a3=970;
  a4=4.05;
  a5=1.8;
  a6=149280;
  a7=5.04;
  a8=130;
  a9=998;
  a10=4.18;
  a11=0.7;
  a12=0.25;
  a13=370;
  a14=0.88;
  a15=285;
  a16=0.05;
  a17=0.0345;
  a18=385.877;
  a19=361.51753;
   

  x1cs = a17/a14;
  x2cs = (a3*a4*a18)/(a14*a6);
  x3cs = (a9*a10*a19)/(a14*a6);	
  u1cs = a16/a16;
  f1cs = ((a5/a12)*x1cs*a1)*(exp((-a2*a3*a4)/(a14*a6*x2cs)));
  f2cs = ((a7*a8)/(a12*a3*a4))*x2cs; 
  f3cs = ((a7*a8)/(a12*a9*a10))*x3cs;
  f4cs = ((a7*a8*a5)/(a12*a11*a3*a4))*x2cs;
  f5cs = ((a5*a7*a8)/(a12*a11*a9*a10))*x3cs;
  f6cs = ((a5*a16)/(a12*a11))*u1cs*x3cs;
  f7cs = ((a9*a10*a5*a16*a15)/(a14*a6*a12*a11))*u1cs; 

  x1c = x(1) + x1cs;
  x2c = x(2) + x2cs;
  x3c = x(3) + x3cs;
  u1c = u(1) + u1cs;
  f1c = ((a5/a12)*x1c*a1)*(exp((-a2*a3*a4)/(a14*a6*x2c)));
  f2c = ((a7*a8)/(a12*a3*a4))*x2c; 
  f3c = ((a7*a8)/(a12*a9*a10))*x3c;
  f4c = ((a7*a8*a5)/(a12*a11*a3*a4))*x2c;
  f5c = ((a5*a7*a8)/(a12*a11*a9*a10))*x3c;
  f6c = ((a5*a16)/(a12*a11))*u1c*x3c;
  f7c = ((a9*a10*a5*a16*a15)/(a14*a6*a12*a11))*u1c;
  f1 = f1c - f1cs;
  f2 = f2c - f2cs;
  f3 = f3c - f3cs;
  f4 = f4c - f4cs;
  f5 = f5c - f5cs;
  f6 = f6c - f6cs;
  f7 = f7c - f7cs;

% Right hand sides of ODEs;	
  sys(1)=-x(1) - f1;
  sys(2)=-x(2) - f2 + f3 + f1;	
  sys(3)=f4 - f5 - f6 + f7;
  
function sys = mdlOutputs(t,x,u)
  a3=970;
  a4=4.05;
  a6=149280;
  a9=998;
  a10=4.18;
  a14=0.88;
   
  sys(1)=x(1);
  sys(4)=a14*x(1);
  sys(2)=x(2);
  sys(5)=(a14*a6*x(2))/(a3*a4);
  sys(3)=x(3);
  sys(6)=(a14*a6*x(3))/(a9*a10);



