s-function error

Hi,

Please, I need some help. While executing the code below in my simulink model, I receive an error message and, I would like to know how to solve it.

The error message is as follows:

Output returned by S-function 'dynamic_model' in 'model_4WID/S-Function' during flag=3 call must be a real vector of length 18

the code :

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

switch flag,

  case 0,
    [sys,x0,str,ts]=mdlInitializeSizes;

  case 1,
    sys=mdlDerivatives(t,x,u);

 
  case 3,
    sys=mdlOutputs(t,x,u);

 
  case { 2, 4, 9 },
    sys = [];

 
  otherwise
    DAStudio.error('Simulink:blocks:unhandledFlag', num2str(flag));

end

function [sys,x0,str,ts]=mdlInitializeSizes

sizes = simsizes;
sizes.NumContStates  = 5;
sizes.NumDiscStates  = 0;
sizes.NumOutputs     = 18;
sizes.NumInputs      = 20;
sizes.DirFeedthrough = 0;
sizes.NumSampleTimes = 1;

sys = simsizes(sizes);
x0  = [20 0 0 0 -5]; % [vx vy r x y] vecteur d'état
str = [];
ts  = [0 0];


function sys=mdlDerivatives(~,x,u)
m=1298.9; Iz=1627; lf=1; lr=1.454; br=1.436; bf=br;

r=x(3);phi=atan2(x(2),x(1));

M=[m 0 0 0 0;0 m 0 0 0;0 0 Iz 0 0;0 0 0 1 0;0 0 0 0 1]; %matrice d'inertie
C=[0 m*r 0 0 0;-m*r 0 0 0 0;0 0 0 0 0;0 0 0 cos(phi) -sin(phi);0 0 0 sin(phi) cos(phi)];
F=[u(9)*cos(u(1))-u(13)*sin(u(1))+u(10)*cos(u(2))-u(14)*sin(u(2))+u(11)*cos(u(3))-u(15)*sin(u(3))+u(12)*cos(u(4))-u(16)*sin(u(4));u(13)*cos(u(1))+u(9)*sin(u(1))+u(14)*cos(u(2))+u(10)*sin(u(2))+u(15)*cos(u(3))+u(11)*sin(u(3))+u(16)*cos(u(4))+u(12)*sin(u(4));u(9)*(lf*sin(u(1))+0.5*bf*cos(u(1)))+u(10)*(lr*sin(u(2))-0.5*bf*cos(u(2)))+u(11)*(-lr*sin(u(3))+0.5*br*cos(u(3)))+u(12)*(-lr*sin(u(4))-0.5*br*cos(u(4)))+u(13)*(lf*sin(u(1))-0.5*bf*cos(u(1)))+u(14)*(lf*sin(u(2))+0.5*bf*cos(u(2)))+u(15)*(-lr*sin(u(3))-0.5*br*cos(u(3)))+u(16)*(-lr*sin(u(4))+0.5*br*cos(u(4)));0;0];

de= M\(F-C*x);

sys(1)=de(1);
sys(2)=de(2);
sys(3)=de(3);
sys(4)=de(4);
sys(5)=de(5);



function sys=mdlOutputs(~,x,u)
                                               
   Cz=0.001; %vertical deflection rate of the tyre  
   Cs=50000;
   epsilon=0.015;
   I=2.1;
   Iz=1627;
   R=0.35;
   Ca=30000;
   m=1298.9;
   lf=1; 
   lr=1.454;
   br=1.436;
   bf=br;
   mu=0.9;
   g=9.81;
   cons=m/(lr+lf);
   h=0.5;
   
   r=x(3);phi=atan2(x(2),x(1));

M=[m 0 0 0 0;0 m 0 0 0;0 0 Iz 0 0;0 0 0 1 0;0 0 0 0 1]; %matrice d'inertie
C=[0 m*r 0 0 0;-m*r 0 0 0 0;0 0 0 0 0;0 0 0 cos(phi) -sin(phi);0 0 0 sin(phi) cos(phi)];
F=[u(9)*cos(u(1))-u(13)*sin(u(1))+u(10)*cos(u(2))-u(14)*sin(u(2))+u(11)*cos(u(3))-u(15)*sin(u(3))+u(12)*cos(u(4))-u(16)*sin(u(4));u(13)*cos(u(1))+u(9)*sin(u(1))+u(14)*cos(u(2))+u(10)*sin(u(2))+u(15)*cos(u(3))+u(11)*sin(u(3))+u(16)*cos(u(4))+u(12)*sin(u(4));u(9)*(lf*sin(u(1))+0.5*bf*cos(u(1)))+u(10)*(lr*sin(u(2))-0.5*bf*cos(u(2)))+u(11)*(-lr*sin(u(3))+0.5*br*cos(u(3)))+u(12)*(-lr*sin(u(4))-0.5*br*cos(u(4)))+u(13)*(lf*sin(u(1))-0.5*bf*cos(u(1)))+u(14)*(lf*sin(u(2))+0.5*bf*cos(u(2)))+u(15)*(-lr*sin(u(3))-0.5*br*cos(u(3)))+u(16)*(-lr*sin(u(4))+0.5*br*cos(u(4)));0;0];


   
   u1=(x(2)+lf*x(3))/(x(1)+0.5*bf*x(3));
   u2=(x(2)+lf*x(3))/(x(1)-0.5*bf*x(3));
   u3=(x(2)-lr*x(3))/(x(1)+0.5*br*x(3));    %% les ui représentent les rapports permettant de calculer les angles de dérives de chaque roue%%
   u4=(x(2)-lr*x(3))/(x(1)-0.5*br*x(3));
   U=[u1 u2 u3 u4];
   
   
                          %------charge verticale:----------
                          
       de= M\(F-C*x);                            
       dde= M\(-C*de); %dérivée de de!
       
       Fz1= cons*(0.5*g*lr-0.5*dde(1)*h-lr*h*dde(2)/bf);
       Fz2= cons*(0.5*g*lr-0.5*dde(1)*h+lr*h*dde(2)/bf);
       Fz3= cons*(0.5*g*lr+0.5*dde(1)*h-lr*h*dde(2)/bf); 
       Fz4= cons*(0.5*g*lr+0.5*dde(1)*h+lr*h*dde(2)/bf);                        
       Fz=[Fz1 Fz2 Fz3 Fz4]';
  
            %%% velocity component in the wheel plane : is the longitunal velocity
       v1=(x(1)+0.5*bf*x(3))*cos(u(1))+(x(2)+lf*x(3))*sin(u(1));
       v2=(x(1)-0.5*bf*x(3))*cos(u(2))+(x(2)+lf*x(3))*sin(u(2));
       v3=(x(1)+0.5*br*x(3))*cos(u(3))+(x(2)-lr*x(3))*sin(u(3));
       v4=(x(1)-0.5*br*x(3))*cos(u(4))+(x(2)-lr*x(3))*sin(u(4));
       V=[v1 v2 v3 v4]';
       
   lamda=zeros(4,1);f=length(lamda);
   omega=zeros(4,1);
   alph=zeros(4,1);
   Re=zeros(4,1);
   S=zeros(4,1);
   dz=zeros(4,1);
   Fs=zeros(4,1);
   Ft=zeros(4,1);
  for i=1:4
      dz(i)=-Cz*Fz(i)+ 0.33*R;
      Re(i)=R-dz(i)./3;
      omega(i)=(-R*u(i+16)+u(i+4))/I;
      S(i)=1+omega(i)*Re(i)/V(i);
      alph(i)=atan(U(i))-u(i);
      lamda(i)= mu*Fz(i)*(1-S(i))*(1-epsilon*V(i)*sqrt((S(i))^2 + (tan(alph(i)))^2))/(2*sqrt((Cs^2)*(S(i))^2 + (Ca^2)*(tan(alph(i))^2)));
      if lamda(i)<1
         f(i)=lamda(i)*(2-lamda(i));
         Fs(i)=Ca*f(i)*tan(alph(i))./(1-S(i));
         Ft(i)=Cs*f(i)*S(i)./(1-S(i));
      elseif lamda(i)>1
          f(i)=1;
          Fs(i)=Ca*f(i)*tan(alph(i))/(1-S(i));
          Ft(i)=Cs*f(i)*S(i)/(1-S(i));
      end
  end
   
       
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%      
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Sorties %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%      
       
sys(1)= x(1);
sys(2)= x(2);
sys(3)= x(3);
sys(4)= x(4);
sys(5)= x(5);
sys(6)=Ft(1);
sys(7)=Ft(2);
sys(8)=Ft(3);
sys(9)=Ft(4);
sys(10)=Fs(1);
sys(11)=Fs(2);
sys(12)=Fs(3);
sys(13)=Fs(4);
sys(14)= phi;
sys(15)=Fz(1);
sys(16)=Fz(2);
sys(17)=Fz(3);
sys(18)=Fz(4);