Matlab Code into Simulink Help

I want to put the following code into Simulink. I tried MATLAB Function block but it gives me many errors. May be a S-Function block will help. I want to have a block with 4 inputs (J, P_D, BARatio, NosBlades) and to have 3 outputs (TT, TQ, TE). Here is the original code

J=0.4;

P_D=0.875; 

BARatio=0.55;  

NosBlades=4;    

%%%%%%%%%%

C2(2)=0.208;
C2(3)=0.241;
C2(4)=0.263;
C2(5)=0.276;
C2(6)=0.279;
C2(7)=0.269;
C2(8)=0.241;
C2(9)=0.184;

%%%%%%%%%

PA(2)=0.888;
PA(3)=1.008;
PA(4)=1.055;
PA(5)=1.060;
PA(6)=1.039;
PA(7)=1.000;
PA(8)=0.948;
PA(9)=0.888;

%%%%%%%%

VV=0;
LA=0.97;
RT=0.045;
LM=0.8*11.3;
LI=0.8*12.67;

%%%%%%%%%

for i=2:9

      C1(i)=C2(i)*BARatio*4.0/(NosBlades*0.5);
      CP=NosBlades*C1(i);                         %Nc/D
      X(i)=i/10;                                  %x/R
      AI(i)=0;                                    %Assume initial Angle of Attack=0
      PR(i)=P_D*PA(i);                            %Pitch/diameter ratio at blade element    
    TD=RT-(RT*0.935)*X(i);
    TC(i)=TD/C1(i);    

      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      counter=1;
      SubConverged=0;
      while and(counter <2000, SubConverged==0)
          counter=counter+1;

          HA=atan(PR(i)/(pi*X(i)))/pi*180-AI(i);  
          H1=tan(HA/180*pi);                      %tan phi
          AA=J/(pi*X(i));                         %tan psi (undisturbed flow angle?)
          EI=AA/H1;                               %Ideal Efficiency
          EA(i)=0.9*EI;                           %Local Efficiency Estimate
          AE=0.0;                                 %tan(gamma) =atan(Cd(i)/Cl(i)) assume =0.0
          KF=X(i)*H1;                             %Lambda parameter for Goldstein correction

          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

          SF=NosBlades/(2.0*KF)-0.5;
          F1=cosh(SF);
          F2=cosh(SF*X(i));
          F3=F2/F1;
          F4=acos(F3);
          K=2.0*F4/pi;

          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
          GK(i)=K;                                

          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

          counter2=0;
          SubSubConverged=0;
          while and(counter2<2000, SubSubConverged==0)
              counter2=counter2+1;
              FA(i)=(1-EI)/(EI+AA*AA/EA(i));                      %axial flow factor
              KT(i)=pi*J^2*X(i)*K*FA(i)*(1+FA(i));                %local dKT/dx
              AE=tan(AE/180*pi);                                  %tan(gamma)
              FT(i)=1-EI*(1+FA(i));                               %a' tangential flow factor

              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

              CL(i)=KT(i)*cos(HA/180*pi);
              CL(i)=CL(i)*4/(pi^2)/(X(i)^2);
              CL(i)=CL(i)/((1-FT(i))^2*(1-H1*AE)*CP);
              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

              F6=0.0107+(AI(i)+1.0)*(-0.0015+AI(i)*(0.0015+0.000965*(AI(i)-1.0)));
              F7=-0.03833+(AI(i)+1.0)*(0.0133+AI(i)*(-0.015-0.01166*(AI(i)-1.0)));
              F8=0.8193+(AI(i)+1.0)*(-0.0138+AI(i)*(0.0903+0.079*(AI(i)-1.0)));
              F9=-3.076+(AI(i)+1.0)*(-0.0728+AI(i)*(-0.3162-0.2437*(AI(i)-1.0)));
              CD(i)=F6+TC(i)*(F7+(TC(i)-0.06)*(F8+F9*(TC(i)-0.12)));
              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

              AG=CD(i)/CL(i);                                     
              AE=atan(AG)/pi*180;                                 %gamma
              EZ=AA/tan((HA+AE)/180*pi);                          %

              if (abs(EZ-EA(i)) < 0.001)                          

                  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

                  K2(2)=1.0+2.857*(BARatio-0.4)*(BARatio-0.6)*(BARatio-0.9);
                  K2(3)=1.19+(BARatio-0.4)*(BARatio-0.6)*(0.267+(BARatio-0.9)*0.1665);
                K2(4)=1.3+(BARatio-0.4)*(0.1+(BARatio-0.6)*(1+(BARatio-0.9)*0.43));
                K2(5)=1.54+(BARatio-0.4)*(0.15+(BARatio-0.6)*(1.1+0.286*(BARatio-0.9)));
                K2(6)=1.67+(BARatio-0.4)*(0.55+(BARatio-0.6)*(0.1667+(BARatio-0.9)*2.9465));
                K2(7)=1.8+(BARatio-0.4)*(0.75+(BARatio-0.6)*(0.3+(BARatio-0.9)*2.835));
                K2(8)=1.8+(BARatio-0.4)*(1.0+(BARatio-0.6)*(1.333+(BARatio-0.9)*1.905));
                K2(9)=1.75+(BARatio-0.4)*(1.25+(BARatio-0.6)*(1.5+(BARatio-0.9)*3.55));

              U1=-0.65*KF*KF+1.1*KF+0.664;

              U2=(0.85+(KF-0.3)*(-4.0+(KF-0.4)*(15.42-47.95*(KF-0.5))));

              U2=-0.09+(KF-0.2)*U2;

              U3=(1.375+(KF-0.3)*(-3.75+(KF-0.4)*(20.85-75.7875*(KF-0.5))));

              U3=-0.2+(KF-0.2)*U3;

              K1=U1+(BARatio-0.4)*(U2+U3*(BARatio-0.8));

              CC(i)=K1*K2(i);

                  if VV==1

                      MC(i)=MT(i)*TC(i);
                      AC=(CC(i)*CL(i)/LI-MC(i))*LM/0.1097+(MC(i)*LI/LA);
                  elseif VV==5
                      MC(i)=0.5*TC(i);
                      AC=(CC(i)*CL(i)/LI-MC(i))*LM/0.1097+(CL(i)/LA);
                      MT(i)=MC(i)/TC(i);
                  else
                      AC=CL(i)/LA;
                      MC(i)=CL(i)/LI;
                      MC(i)=MC(i)*CC(i);

                      if (MC(i)<0.5*TC(i))
                          MT(i)=MC(i)/TC(i);
                      else
                          VV=5;
                          MC(i)=0.5*TC(i);
                          AC=(CC(i)*CL(i)/LI-MC(i))*LM/0.1097+(CL(i)/LA);
                          MT(i)=MC(i)/TC(i);
                      end              

                  end

                 SubSubConverged=1;                                %Set subsub loop converged flag to 1

              else
                  EA(i)=EZ;                                        %update local efficiency estimate
              end          
          end %subsubloop
          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

          if (abs((AC-AI(i))/(AC))<0.1) 
                             %subloop converged difference between angle of attack assumed AI and calculated angle of attack AC <0.1*AC
              KQ(i)=4.935*J*X(i)*X(i)*X(i)*K*FT(i)*(1+FA(i));      %calculate dkq/dx
              SubConverged=1;                                      %set subloop converged flag to 1
          else
              AI(i)=(AC+AI(i))/2;                                  %New Guess of Angle of Attack AI(i)
          end
      end %subloop

      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  end     %End Main Loop

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

X(10)=1.0;

KT(10)=0.0;

KQ(10)=0.0;

FA(10)=0.0;

FT(10)=0.0;

GK(10)=0.0;

MT(10)=MT(9);

PR(10)=PR(9);

AI(10)=AI(9);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Numerical Integration of Total Values

TT=(KT(2)+2.0*(KT(4)+KT(6)+KT(8))+4.0*(KT(3)+KT(5)+KT(7)+KT(9)));   

TT=0.1*TT/3.0;                                                      

TQ=(KQ(2)+2.0*(KQ(4)+KQ(6)+KQ(8))+4.0*(KQ(3)+KQ(5)+KQ(7)+KQ(9)));   

TQ=0.1*TQ/3.0;                                                      

TE=J*TT/(2.0*3.1416*TQ);