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

global c p i
switch flag
    
    case 0  %Initialize
      sys = [3+2*i.n,      % number of continuous states
            0,       % number of discrete states
            3+2*i.n,       % number of outputs
            5,       % number of inputs
            0,       % reserved must be zero
            0,       % direct feedthrough flag
            1];      % number of sample times
        
        x0 = [i.N_1_s ones(1,i.n/2+1)*i.T_1 ones(1,i.n/2+1)*32 ...
            ones(1,i.n/2)*c.T_cin ones(1,i.n/2)*c.T_cin]';
        str = [];
        ts  = [0 0];   % sample time: [period, offset]
        
    case 1 %Derivatives
        
        % State equations:
             %Solids
        dxdt(1) = u(1)-u(2);  %Mass balance top section
        dxsdt(1) = (1/x(1))*(u(3)-x(2))*u(1); %Energy wrong top section 
        s = [];  %Stage counter
        
        xl(1)= 0.0;         %Mock state for retaining index
        xl(i.n/2+2)=c.T_cin;%matching the solid phase
        for j = 2:1:i.n+3 %Extracting the state variables for solid EB
               xs(j-1) = x(j);    
               if  j<i.n/2+3 & j>2
                    xl(j-1) = x(j+i.n+1);   %Counter current
               elseif j>i.n/2+3             %Skip middle section
                   xl(j-1) = x(j+i.n);      %Co-current
               end
        end
                    % EB Solid - Counter current Hx
                    % i.n = 10 -> Section 2-6
        for s = 2:1:i.n/2+1
            if x(s) > 34
                p.a_NPK1 =p.a_NPKtop;
            elseif x(s) > 33
                 p.a_NPK1 =p.a_NPKphase;
            else
                p.a_NPK1 =p.a_NPKbot;
            end
              dxsdt(s) = (1/i.N_i_s)*(xs(s-1)-xs(s))*u(2)...
                     -(p.U*p.Ai/(i.N_i_s*p.a_NPK1))*(xs(s)-xl(s));
                    
                    % Transfer back to program state vector
                    dxdt(s+1) = dxsdt(s);
        end
        
                    % EB Solid - Middle section
                    % i.n = 10 -> Section 7
        dxdt(i.n/2+3)= (1/i.N_7_s)*(xs(i.n/2+1)-xs(i.n/2+2))*u(2); 
        
                    % EB Solid - Co-current Hx
                    % i.n = 10 -> Section 8-12
        s = [];
        for s = (i.n/2+3):1:i.n+2 
         if x(s) > 34
                p.a_NPK1 =p.a_NPKtop;
            elseif x(s) > 33
                 p.a_NPK1 =p.a_NPKphase;
            else
                p.a_NPK1 =p.a_NPKbot;
            end
             dxsdt(s) = (1/i.N_i_s)*(xs(s-1)-xs(s))*u(2)...
            -(p.U*p.Ai/(i.N_i_s*p.a_NPK1))*(xs(s)-xl(s));
             
            dxdt(s+1) = dxsdt(s);
        end
                    % EB Liquid - Co current HX
                    % i.n = 10 -> Section 8-12
        for s = (i.n/2+3):1:i.n+2
                dxldt(s) = (1/i.N_i_l)*(xl(s-1)-xl(s))*u(4)+...
                (p.U*p.Ai/(i.N_i_l*p.a_H2O1))*(xs(s)-xl(s));
                dxdt(s+1+i.n)=  dxldt(s);
        end
                    % EB Liquid - Counter current HX
                    % i.n = 10 -> Section 6
        s = i.n/2+1;
        dxldt(s) = (1/i.N_i_l)*(xl(i.n+2)-xl(s))*u(4)...
                 +(p.U*p.Ai/(i.N_i_l*p.a_H2O1))*(xs(s)-xl(s));
                    dxdt(s+2+i.n)=  dxldt(s);
              
                    % EB Liquid - Counter current HX
                    % i.n = 10 -> Section 2-5
        for s = i.n/2:-1:2 
                dxldt(s) = (1/i.N_i_l)*(xl(s+1)-xl(s))*u(4)...
                 +(p.U*p.Ai/(i.N_i_l*p.a_H2O1))*(xs(s)-xl(s));
                dxdt(s+2+i.n)= dxldt(s);
        end
%                     s = 6
%                     disp('start')
%                     disp(dxldt(s))
%                     disp(xl(12))
%                     disp(xl(6))
%                     disp((1/(i.N_i_s)/5)*(xl(i.n+2)-xl(s))*u(4))
%                     disp((p.U*p.Ai/((i.N_i_s)/5*p.a_H2O1))*(xs(s)-xl(s)))
%         

        
%         disp(dxldt)
%         disp(dxsdt)
%      disp('xl')
%     disp(xl')
%      disp(dxldt(6))
%         disp('xs')
%     disp(xs')
        dxdt = dxdt(:);
        sys = dxdt;
        
    case 2 % Discrete state udate
        
        sys = []; %There is none
        
    case 3 % Outputs
        
        j= [];
        for j=1:1:i.n*2+3
            sys(j) = [x(j)];
        end
        sys = sys';
        
    case 9 % Terminate
        
        sys = [];
        
    otherwise
        error(['unhandled flag = ',num2str(flag)]);
        
end
        