Modelling a discrete system using level 1 matlab s function.

Hello everyone, I have modeled a gas pipeline to find out its output pressure and input gas flow using the code given below. The system has two o/p and two i/p and it is solved numerically. Now i want to model the gas duct using level 1 matlab s function block and am confused how to put the code below into the s function. Is there any method of handling numerical methods in a s function.??? Please help!! The m file is attached herewith. Thank you!

g=9.81; Qo=50; f=.003; d=0.73; D=0.6; S=3.14*D^2/4; L= 10^5;R=392; T= 278; P0 = 50e5; dx=5000; dt=.25/60; 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 disp(Z) c=sqrt(Z*R*T); disp(c) k1 = (c*d)/S; k2 = (2*f/D)*k1^2*(dx/2); k3 = (g*sin(theta)*dx)/(2*c^2);

T= 2*24; %No of time steps N=(T/dt)+1; I=(L/dx)+1;%No of distance steps

P_A= zeros(N,I); Q_A=zeros(N,I); P_B= zeros(N-1,I-1); Q_B=zeros(N-1,I-1);

t=0:dt:T; xx=0:dx:L; P_A(:,1)= 50e5; Q_A(1,:)=50; Q_A(:,I)=interp1(Time,Demand,t); % input

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(1,i)=init_press(xx(i)); end

for n=1:N % stating to solve the o/p for each time step n1=n; if n1~=N i=1; % calculate layer B values for i1=1:I-1 d1=(1+0.5*k3)*P_A(n,i+1) - k1*Q_A(n,i+1) + (0.5*k2)*(Q_A(n,i+1)^2/P_A(n,i+1)); d2=(0.5*k3-1)*P_A(n,i) - k1*Q_A(n,i) + (0.5*k2)*(Q_A(n,i)^2/P_A(n,i)); P_B(n1,i1)=0.5*(d1-d2);

            a=0.5*k2; b=k1*P_B(n1,i1); 
            c=d2*P_B(n1,i1) + (1+0.5*k3)*P_B(n1,i1)^2;

            Q_B(n1,i1)= (-b + sqrt(b^2 - 4*a*c))/(2*a);
            i=i+1;
        end
       % produce next layer A values using Layer B values
         a=0.5*k2; b=k1*P_A(n+1,1); 
         d1=(1+0.5*k3)*P_B(n1,1) - k1*Q_B(n1,1) + (0.5*k2)*(Q_B(n1,1)^2/P_B(n1,1));
            c=d1*P_A(n+1,1) + (0.5*k3-1)*P_A(n+1,1)^2;
        Q_A(n+1,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(n1,i1+1) - k1*Q_B(n1,i1+1) + (0.5*k2)*(Q_B(n1,i1+1)^2/P_B(n1,i1+1));
            d2=(0.5*k3-1)*P_B(n1,i1) - k1*Q_B(n1,i1) + (0.5*k2)*(Q_B(n1,i1)^2/P_B(n1,i1));
            P_A(n+1,i)=0.5*(d1-d2);

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

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

            i1=i1+1;
        end

        i=I;
        a=0.5*k3+1; 
        d2=(0.5*k3-1)*P_B(n1,I-1) - k1*Q_B(n1,I-1) + (0.5*k2)*(Q_B(n1,I-1)^2/P_B(n1,I-1)); %P_B(n1,I-1) & Q_B(n1,I-1) are the states
        b=k1*Q_A(n+1,I) + d2;
        c=0.5*k2*Q_A(n+1,I)^2; % Q_A(n+1,I) - current input
        P_A(n+1,I) = (-b + sqrt(b^2 - 4*a*c))/(2*a); % Output 2
    end
end
P_out= P_A(:,I)./10^5; %the outputs from s function block 
Q_out= Q_A(:,1);