Why simulation is taking too long to respond?

Hello everyone..I am quite new to writing matlab s functions.Still i have written a matlab level-2 S function to model a physical system and eventually it has become large cause there are too many calculations. Ok..now while i am simulating the model including the s func block...the simulation is taking too long to respond and its just showing the word 'Initializing' in the bottom left corner and after sometime matlab stops responding. Why is this happening? Would u please check my code and tell whats wrong? The code is attached. Thank you!!

function moc_sfunc1(block)

setup(block);

%endfunction

function setup(block)

%%Register dialog parameter:
block.NumDialogPrms = 0;

%%Register number of input and output ports
block.NumInputPorts  = 10;
block.NumOutputPorts = 2;

%%Setup functional port properties to dynamically inherited.
block.SetPreCompInpPortInfoToDynamic;
block.SetPreCompOutPortInfoToDynamic;

block.InputPort(1).Complexity   = 'Real'; 
block.InputPort(1).DataTypeId   = 0; %double
block.InputPort(1).SamplingMode = 'Sample';
block.InputPort(1).Dimensions   = 1;

block.InputPort(2).Complexity   = 'Real';
block.InputPort(2).DataTypeId   = 0;
block.InputPort(2).SamplingMode = 'Sample';
block.InputPort(2).Dimensions   = 1;

block.InputPort(3).Complexity   = 'Real';
block.InputPort(3).DataTypeId   = 0;
block.InputPort(3).Dimensions   = 1;
block.InputPort(3).SamplingMode = 'Sample';

block.InputPort(4).Complexity   = 'Real';
block.InputPort(4).DataTypeId   = 0;
block.InputPort(4).Dimensions   = 1;
block.InputPort(4).SamplingMode = 'Sample';

block.InputPort(5).Complexity   = 'Real';
block.InputPort(5).DataTypeId   = 0;
block.InputPort(5).Dimensions   = 1;
block.InputPort(5).SamplingMode = 'Sample';

block.InputPort(6).Complexity   = 'Real';
block.InputPort(6).DataTypeId   = 0;
block.InputPort(6).Dimensions   = 1;
block.InputPort(6).SamplingMode = 'Sample';

block.InputPort(7).Complexity   = 'Real';
block.InputPort(7).DataTypeId   = 0;
block.InputPort(7).Dimensions   = 1;
block.InputPort(7).SamplingMode = 'Sample';

block.InputPort(8).Complexity   = 'Real';
block.InputPort(8).DataTypeId   = 0;
block.InputPort(8).Dimensions   = 1;
block.InputPort(8).SamplingMode = 'Sample';

block.InputPort(9).Complexity   = 'Real';
block.InputPort(9).DataTypeId   = 0;
block.InputPort(9).Dimensions   = 1;
block.InputPort(9).SamplingMode = 'Sample';

block.InputPort(10).Complexity   = 'Real';
block.InputPort(10).DataTypeId   = 0;
block.InputPort(10).Dimensions   = 1;
block.InputPort(10).SamplingMode = 'Sample';

block.OutputPort(1).Complexity   = 'Real';
block.OutputPort(1).DataTypeId   = 0;
block.OutputPort(1).SamplingMode = 'Sample';
block.OutputPort(1).Dimensions   = 1;

block.OutputPort(2).Complexity   = 'Real';
block.OutputPort(2).DataTypeId   = 0;
block.OutputPort(2).SamplingMode = 'Sample';
block.OutputPort(2).Dimensions   = 1;

% Register the sample times.
block.SampleTimes = [-1 0];

% -----------------------------------------------------------------
% Options
% -----------------------------------------------------------------
% Specify if Accelerator should use TLC or call back to the 
% MATLAB file
block.SetAccelRunOnTLC(false);

%%Set the block simStateCompliance to default (i.e., same as a built-in block)
block.SimStateCompliance = 'DefaultSimState';

%%Register methods %%

block.RegBlockMethod('PostPropagationSetup',    @DoPostPropSetup);
block.RegBlockMethod('InitializeConditions',    @InitializeConditions);
block.RegBlockMethod('Start',                   @Start);  
block.RegBlockMethod('Outputs',                 @Outputs);

%endfunction

function DoPostPropSetup(block)

%%Setup Dwork  
I = 21;                    %%no. of length steps  
block.NumDworks = 7;
block.Dwork(1).Name = 'Q_A'; %Q_A(1),...Q_A(I)
block.Dwork(1).Dimensions      = I;
block.Dwork(1).DatatypeID      = 0;
block.Dwork(1).Complexity      = 'Real';
block.Dwork(1).UsedAsDiscState = true;

block.Dwork(2).Name = 'P_A'; %%P_A(1),....P_A(I)
block.Dwork(2).Dimensions      = I;
block.Dwork(2).DatatypeID      = 0;
block.Dwork(2).Complexity      = 'Real';
block.Dwork(2).UsedAsDiscState = true;

block.Dwork(3).Name = 'Q_B'; %Q_B(1),...Q_B(I-1)
block.Dwork(3).Dimensions      = I-1;
block.Dwork(3).DatatypeID      = 0;
block.Dwork(3).Complexity      = 'Real';
block.Dwork(3).UsedAsDiscState = true;

block.Dwork(4).Name = 'P_B';  %P_B(1),....P_B(I-1)
block.Dwork(4).Dimensions      = I-1;
block.Dwork(4).DatatypeID      = 0;
block.Dwork(4).Complexity      = 'Real';
block.Dwork(4).UsedAsDiscState = true;

block.Dwork(5).Name = 'k1';  
block.Dwork(5).Dimensions      = 1;
block.Dwork(5).DatatypeID      = 0;
block.Dwork(5).Complexity      = 'Real';
block.Dwork(5).UsedAsDiscState = true;

block.Dwork(6).Name = 'k2'; 
block.Dwork(6).Dimensions      = 1;
block.Dwork(6).DatatypeID      = 0;
block.Dwork(6).Complexity      = 'Real';
block.Dwork(6).UsedAsDiscState = true;

block.Dwork(7).Name = 'k3'; 
block.Dwork(7).Dimensions      = 1;
block.Dwork(7).DatatypeID      = 0;
block.Dwork(7).Complexity      = 'Real';
block.Dwork(4).UsedAsDiscState = true;

%%Register all tunable parameters as runtime parameters.

%endfunction

function Start(block)

%%Initialize Dwork 
I=21;
block.Dwork(1).Data = zeros(1, I);
block.Dwork(2).Data = zeros(1, I);
block.Dwork(3).Data = zeros(1, I-1);
block.Dwork(4).Data = zeros(1, I-1);

%endfunction

function InitializeConditions(block)

g=9.81;

Qo=block.InputPort(8).Data;

f=block.InputPort(5).Data;

d=block.InputPort(3).Data;

D=block.InputPort(4).Data;

S=3.14*D^2/4;

L= block.InputPort(7).Data;

R=392;

T=block.InputPort(10).Data;

P0 =block.InputPort(9).Data;

H=block.InputPort(6).Data;

dx=5000;

theta=atan(H/L); % all the parameter values

Z0=.8;

CF=0;

while CF==0

c=sqrt(Z0*R*T); % a loop to find out the actual value of compressibility

if H~=0

sigma=(2*g*sin(theta))/c^2;

lamda=(2*Qo^2*f*c^4*d^2)/(D*S^2*g*sin(theta));

PL = sqrt((P0^2-lamda*(exp(sigma*L)-1))*exp(-sigma*L));

else phy=(f/D)*(2*c*d*Qo/S)^2;

    PL=sqrt(P0^2-phy*L);
end

Pav= 0.5*(P0+PL);

Z= 1-(Pav/390e5);

if abs(Z-Z0)<10^-6

    CF=1;

end

Z0=Z;

end

c=sqrt(Z*R*T);

k1 = (c*d)/S;

k2 = (2*f/D)*k1^2*(dx/2);

k3 = (g*sin(theta)*dx)/(2*c^2);

 Q_A = block.Dwork(1).Data;
 P_A = block.Dwork(2).Data;

 Q_A(1:I) = Qo;
 xx=0:dx:L;
 if H~=0
init_press=@(x)(sqrt((P0^2-lamda*(exp(sigma*x)-1))*exp(-sigma*x)));
else
    init_press=@(x)(sqrt(P0^2-phy*x));
end
for i=1:I
P_A(i)=init_press(xx(i)); %initial pressure condition
end
 block.Dwork(1).Data = Q_A;
 block.Dwork(2).Data = P_A;
 block.Dwork(5).Data = k1;
 block.Dwork(6).Data = k2;
 block.Dwork(7).Data = k3;

function Outputs(block)

k1 = block.Dwork(5).Data; k2 = block.Dwork(6).Data; k3 = block.Dwork(7).Data;

 Q_A = block.Dwork(1).Data;
 P_A = block.Dwork(2).Data;
 Q_B = block.Dwork(3).Data;
 P_B = block.Dwork(4).Data;

 P_in = block.InputPort(2).Data;
 Q_dem = block.InputPort(1).Data;

i=1; % calculate layer B values

        for i1=1:I-1

            d1=(1+0.5*k3)*P_A(i+1) - k1*Q_A(i+1) + (0.5*k2)*(Q_A(i+1)^2/P_A(i+1));

            d2=(0.5*k3-1)*P_A(i) - k1*Q_A(i) + (0.5*k2)*(Q_A(i)^2/P_A(i));

            P_B(i1)=0.5*(d1-d2);

            a=0.5*k2; b=k1*P_B(i1); 

            c=d2*P_B(i1) + (1+0.5*k3)*P_B(i1)^2;

            Q_B(i1)= (-b + sqrt(b^2 - 4*a*c))/(2*a);

            i=i+1;

        end
       % produce current layer A values using Layer B values
         a=0.5*k2; b=k1*P_in; 
         d1=(1+0.5*k3)*P_B(1) - k1*Q_B(1) + (0.5*k2)*(Q_B(1)^2/P_B(1));
            c=d1*P_in + (0.5*k3-1)*P_in^2;
        Q_A(1)= (-b + sqrt(b^2 - 4*a*c))/(2*a); %  Output 1
        i1=1;
        for i=2:I-1
            d1=(1+0.5*k3)*P_B(i1+1) - k1*Q_B(i1+1) + (0.5*k2)*(Q_B(i1+1)^2/P_B(i1+1));
            d2=(0.5*k3-1)*P_B(i1) - k1*Q_B(i1) + (0.5*k2)*(Q_B(i1)^2/P_B(i1));
            P_A(i)=0.5*(d1-d2);

            a=0.5*k2; b=k1*P_A(i); 
            c=d2*P_A(i) + (1+0.5*k3)*P_A(i)^2;

            Q_A(i)= (-b + sqrt(b^2 - 4*a*c))/(2*a);

            i1=i1+1;
        end

        a=0.5*k3+1; 
        d2=(0.5*k3-1)*P_B(I-1) - k1*Q_B(I-1) + (0.5*k2)*(Q_B(I-1)^2/P_B(I-1)); 
        b=k1*Q_dem + d2;
        c=0.5*k2*Q_dem^2; % Q_dem - current input
        P_A(I) = (-b + sqrt(b^2 - 4*a*c))/(2*a); % Output 2

block.Dwork(1).Data = Q_A;
block.Dwork(2).Data = P_A;
block.Dwork(3).Data = Q_B;
block.Dwork(4).Data = P_B;

block.OutputPort(1).Data = Q_A(1);
block.OutputPort(2).Data = P_A(I);

%endfunction

[roro6992]

Please tell me if there is any wrong in my s function especially in the initializeconditions method. Thank you

[Geoff Hayes]

I don't have simulink so can't run the above, but if the code is getting stuck while intializing you may want to look at your InitializeConditions function. I noticed that there is a while loop within it - perhaps the code is getting stuck here.

Can you add some fprintfs to write out values for Z,*Z0*, and CF at each iteration?