what is wrong with my code in Matlab function?

Hi, please i dont know what is wrong with my code. It is gving error as follows: ''Error in default port dimensions function of S-function 'space_time_kintetics_model_trial/MATLAB Function'. This function does not fully set the dimensions of output port 2'

I want to model a system such that I can get the total output as well as output from different sections of the segmented volume of the system. It means i need 5 outputs. the characteristic equations for different segment differs. So i used the code as shown below and the attached Matlab FCN file. Kinldy help to improve the code or the likely cause of the error. I want to further apply a controller to the system to form a closed-loop for further analysis. I have applied the matlab file of the system. Thank you

function [dnr,dnr1,dnr2,dnr3,dnr4,dCr1,dCr2,dCr3,dTf,dTc,dX,dI,drho_r1,drho_r2,drho_r3,drho_r4,dh1,dh2,rho] = reactor_space_time_kinetics(nr,Cr1,Cr2,Cr3,Tf,Tc,X,I,h1,h2,Zr1,Zr2,rho_r1,rho_r2,rho_r3,rho_r4)
P0 = 3000;
h=355;%height of the core
d=316;%diameter of the core 
beta=0.0065;
Tcin=291;%coolant inlet temperature 
% w=16704;%coolant flow rate 
D = 0.16;% diffusion constant 
v = 22000;%mean velocity of the thermal neutron
l = 2e-5;% effective prompt neutron lifetime 
S=3.14*d^2/4; %surfcae interface of node
V=(3.14*h*d^2/4)/4;% volume of each node
Sig_x = 2.36e-18;%Microscopic absorption cross section of xenon X
beta1 = 0.0002145;% delayed neutron fraction of group 1 precurssor
beta2 = 0.00225;% delayed neutron fraction of group 2 precurssor
beta3 = 0.00404;% delayed neutron fraction of group 3 precurssor
Lemda_1 = 0.0124;% Radiaoctive decay constant of first group neutron precussor 
Lemda_2 = 0.0369;% Radiaoctive decay constant of second group neutron precussor 
Lemda_3 = 0.632;% Radiaoctive decay constant of third group neutron precussor 
Lemda_x = 2.08e-5;% decay constant of xenon
Lemda_I = 2.88e-5;% decay constant of iodine
gama_x = 0.003;% fractional yield of xenon
gama_I = 0.059;%fractional yield of iodine
Gr = 0.0145; % Total reactivity rod worth
Sum_f = 0.3358;%Macroscopic fission cross section
f_f=0.92;
G=3.2e-11;
h0=355; %total core height
nr0 = 1;
Tf0=30;
Tc0=30;
X0=1;
P=P0*nr;
M=(280*nr0+74);
alfa_f=(nr0-4.24)*1e-5;
mu_f=26.3;
alfa_c=(-4*nr0-17.3)*1e-5;
mu_c=(160*nr0/9+54.022);
Omega=(5*nr0/3+4.93333);
drho_r1=Gr*Zr1/2*(1-sign(h1-h1/4));
drho_r2=Gr*Zr1/2*(1-sign(h1-h1/4));
drho_r3=Gr*Zr2/2*(1-sign(h2-3*h2/4));
drho_r4=Gr*Zr2/2*(1-sign(h2-3*h2/4));
dh1=h0*Zr1;
dh2=h0*Zr2;
rho_r=rho_r1+rho_r2+rho_r3+rho_r4;
rho_fT=alfa_f*(Tf-Tf0)+alfa_c*(Tc-Tc0);
rho_fp=Sig_x*(X-X0)/v*Sum_f;
rho=rho_r+rho_fT+rho_fp;
dnr= (rho-beta)/l*nr+beta1/l*Cr1+beta2/l*Cr2+beta3/l*Cr3;
dCr1 = Lemda_1*nr-Lemda_1*Cr1;
dCr2= Lemda_2*nr-Lemda_2*Cr2;
dCr3= Lemda_3*nr-Lemda_3*Cr3;
dI= gama_I*P/G*V-Lemda_I*I;
dX= gama_x*P/G*V+Lemda_I*I-(Sig_x*P/G*V*Sum_f+Lemda_x)*X;
dTf =1/mu_f*( f_f*P-Omega*(Tf-Tc));
dTc =1/mu_c*(1-f_f)*P+Omega*(Tf-Tc-2*M*(Tc-Tcin));
alfa=D*v*l*S/d*V;
% alfa=sum(alfa);
% rho(i)= rho_r(i)+alfa_f(i)*(Tf1-30)+alfa_c(i)*(Tc1-30)-Sig_x*(X1-1)/v*Sum_f;

 
 for i=1:4 
        for j=1:4 
             if j~=i
alfa(j,i)=D(i)*v*l(i)*S(i)/d(i,j)*V(i);
             else 
alfa(i,i)=sum(alfa(j,i));
             end
dnr(i) = (rho(i)-beta)/l(i)*nr(i)+beta1/l(i)*Cr1(i)+beta2/l(i)*Cr2(i)+beta3/l(i)*Cr3(i)-alfa(i,i)*nr(i)/l(i)+alfa(j,i)*nr(j)/l(i);
             
        end
dCr1(i) = Lemda_1*nr(i)-Lemda_1*Cr1(i);
dCr2(i) = Lemda_2*nr(i)-Lemda_2*Cr2(i);
dCr3(i) = Lemda_3*nr(i)-Lemda_3*Cr3(i);

dI(i) = gama_I*P(i)/G*V(i)-Lemda_I*I(i);
dX(i) = gama_x*P(i)/G*V(i)+Lemda_I*I(i)-(Sig_x*P(i)/G*V(i)*Sum_f+Lemda_x)*X(i);
dTf(i) =1/mu_f(i)*( f_f*P(i)-Omega(i)*(Tf(i)-Tc(i)));
dTc(i) =1/mu_c(i)*(1-f_f)*P(i)+Omega(i)*(Tf(i)-Tc(i)-2*M(i)*(Tc(i)-Tcin(i)));
mu_c(i)=(160*nr0(i)/9+54.022);
Omega(i)=(5*nr0(i)/3+4.93333);
M(i)=(280*nr0(i)+74);
rho_fT(i)=alfa_f(i)*(Tf(i)-Tf0(i))+alfa_c(i)*(Tc(i)-Tc0(i));
alfa_f(i)=(nr0(i)-4.24)*1e-5;
alfa_c(i)=(-4*nr0(i)-17.3)*1e-5;
rho_fp(i)=Sig_x*(X(i)-X0(i))/v*Sum_f;
% rho(i)= rho_r(i)+alfa_f(i)*(Tf1-30)+alfa_c(i)*(Tc1-30)-Sig_x*(X1-1)/v*Sum_f;
rho(i)=rho_r(i)+rho_fT(i)+rho_fp(i);
P(i)=P0(i)*nr(i);
nr=sum(nr(i))/4;
% dnr(i)=dnri;
% dnr(2)=dnr2;
% dnr(3)=dnr3;
% dnr(4)=dnr4;
% nr=nr(1)+nr(2)+nr(3)+nr(4)/4;


 end 

[Stephen23]

Rather than lots of position inout/output arguments, it would be likely simpler and more robust to pass those parameters in one simple scalar structure.

[KaMATLAB]

@Stephen Cobeldick thank you for the insight

Answer 1

[Cris LaPierre]

See this post.