Drum-Boiler model script not running properly

I am trying to run a simulation of a Drum-Boiler modeled by Astrom and Bell (1998). I have used the same code available in their paper which uses a S-function with 4 inputs and 11 outputs. My problem is because I need some variables that are calculated in flag 1 and used in flag 3. So the script is calling flag 3 prior to flag 1 so my variables are empty.

Could anyone help me to solve this?

Attached there is my Simulink file and the article of the model.

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

%sfd302 S-function version of the four state model

%compatible with simulink block diagrams.

global qr qdc qct;

% qr = 0;

% qdc = 0;

% qct = 0;

%plant parameters-------------------------------

mt=300000;mr=100000;md=20000;Cp=650;

Cfw=4.18;Ad=23;Vd=37;Vr=37;Vdc=11;

Vt=Vd+Vr+Vdc;

ke=200;Vsd0=8;beta=0.3;

%properties of steam and water in saturated table

a01=2.7254E6;

a11=-1.8992E4;

a21=-1160.0;

a02=53.1402;

a12=7.673;

a22=0.36;

a03=1.4035E6;

a13=4.9339E4;

a23=-880.0;

a04=691.35;

a14=-18.672;

a24=-0.0603;

a05=310.6;

a15=8.523;

a25=-0.33;

%-------------------------------------------------

switch flag,

case 0, %INITIALIZATIONS

str=[];

sizes = simsizes;

sizes.NumContStates = 4;

sizes.NumDiscStates = 0;

sizes.NumOutputs = 11;

sizes.NumInputs = 4;

sizes.DirFeedthrough = 0;

sizes.NumSampleTimes = 1;

sys = simsizes(sizes);

ts=[0,0];

x0(1)=x_init(1);

x0(2)=x_init(2);

x0(3)=x_init(3);

x0(4)=x_init(4);

case 1, %compute state derivatives

%state initial conditions

Vwt=x(1);

p=x(2);

ar=x(3);

Vsd=x(4);

%controls or inputs

Q=u(1)*1e6; %conversion to watts

qf=u(2);

tf=u(3);

qs=u(4)+4.8*(p-8.5); %correction for steam transducer

% properties of steam and water in saturated state

hs=a01+(a11+a21*(p-10))*(p-10);

dhsdp=a11+2*a21*(p-10);

rs=a02+(a12+a22*(p-10))*(p-10);

drsdp=a12+2*a22*(p-10);

hw=a03+(a13+a23*(p-10))*(p-10);

dhwdp=a13+2*a23*(p-10);

rw=a04+(a14+a24*(p-10))*(p-10);

drwdp=a14+2*a24*(p-10);

ts=a05+(a15+a25*(p-10))*(p-10);

dtsdp=a15+2*a25*(p-10);

%properties of water in subcritical state

hf=(Cfw*tf+p*1e3/rw)*1e3;

%the model equations--------------------

hc=hs-hw;

hr=ar*hs+(1-ar)*hw;

%average steam quality volume rati and partial derivatives

z=ar*(rw-rs)/rs;

av=rw/(rw-rs)*(1-(1/z)*log(1+z));

davdar=rw/rs*((1/z)*log(1+z)-1/(z+1))/z;

zp=(rs*drwdp-rw*drsdp)*ar/rs/rs;

z1=(rw*drsdp-rs*drwdp)/(rw-rs)/(rw-rs);

z2=rw/(rw-rs)*((1/z)*log(1+z)-1/(1+z))/z*zp;

davdp=z1*(1-(1/z)*log(1+z))+z2;

%circulation flow

s1=ke*(rw-rs)*Vr*av;

qdc=sqrt(s1);

%equations for coefficents of derivatives of state variables

Td=600/qs;

Vst=Vt-Vwt;

Vwd=Vwt-Vdc-(1-av)*Vr;

e11=rw-rs;

e12=Vst*drsdp+Vwt*drwdp;

e21=hw*rw-hs*rs;

e22x=-Vt*1e6+mt*Cp*dtsdp;

e22=Vst*(hs*drsdp+rs*dhsdp)+Vwt*(hw*drwdp+rw*dhwdp)+e22x;

e3w=(rw*dhwdp-ar*hc*drwdp)*(1-av)*Vr;

e3x=(rs+(rw-rs)*ar)*hc*Vr*davdp-Vr*1e6+mr*Cp*dtsdp;

e32=((1-ar)*hc*drsdp+rs*dhsdp)*av*Vr+e3w+e3x;

e33=(rs+(rw-rs)*ar)*hc*Vr*davdar;

e42y=(1+beta)*ar*Vr;

e42x=e42y*((1-av)*drwdp+av*drsdp+(rs-rw)*davdp)+Vsd*drsdp;

e42=(rs*Vsd*dhsdp+rw*Vwd*dhwdp-Vsd*1e6+md*Cp*dtsdp)/hc+e42x;

e43=(1+beta)*(ar*Vr*(rs-rw)*davdar);

e44=rs;

e=[e11,e12,0,0

e21,e22,0,0

0,e32,e33,0

0,e42,e43,e44];

%the right hand side of state equations

b4=(hf-hw)*qf/hc+rs*(Vsd0-Vsd)/Td;

b=[qf-qs; Q+qf*hf-qs*hs; Q-qdc*ar*hc; b4];

%solve linear equation for derivatives

dx=e\b;

qx=(rw-rs)*Vr*davdp*dx(2)+(rw-rs)*Vr*davdar*dx(3);

%two important flows for understanding behaviour

qr=qdc-(av*drsdp+(1-av)*drwdp)*Vr*dx(2)+qx;

qctx=rs*Vst*dhsdp+rw*Vwt*dhwdp-Vt*1e6+mt*Cp*dtsdp;

qct=((hw-hf)*qf+qctx*dx(2))/hc;

%step the derivatives, dx is a 4X1 vector

sys = [dx];

%-------------------------------------------------------------

case 3, %Compute outputs

%extract state variables

Vwt=x(1);

p=x(2);

ar=x(3);

Vsd=x(4);

%drum level

rs=a02+(a12+a22*(p-10))*(p-10);

rw=a04+(a14+a24*(p-10))*(p-10);

z=ar*(rw-rs)/rs;

av=rw/(rw-rs)*(1-(1/z)*log(1+z));

Vwd=Vwt-Vdc-(1-av)*Vr;

ls=Vsd/Ad;

lw=Vwd/Ad;

l=lw+ls;

sys=[l, p, lw, ls, qr, av, qdc, ar, Vwt, Vsd, qct];

%--------------------------------------------------------------

case {2,4,9},

sys=[];

otherwise

error(['Unhandled flag = ',num2str(flag)]);

end