Help with formulating conditional logic.

I am writing a stability and control code (for a class assignment) that takes into account 3 different airplane conditions; Climb data, Cruise data, and Approach data. I don't want to change variables for each condition, rather I have captured all of the data for each condition in the code with differing data for each condition.

Ultimately, I would like to output the transfer functions for each condition. I am not very savy with MATLAB and would love some help in finding the best way to formulate this logic.

The code:

clc 
clear all

g = 32.174
theta1 = 0
%Cessna design parameters
S = 174 %Planform Area (ft^2)
c = 4.9 %Chord length (ft)
b = 36.0 %Wing span (ft)

%%%%%%%%%%%%%%%%Flight Conditions%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Climb Flight Condition
altitude = 0
M = 0.120
U_1 = 133.5
q_bar = 21.2
c_frac = 26.2
a_1 = 5.4

%Cruise Flight Condition
h = 5000 %Altitude (ft)
M = 0.201 %Mach Number
U_1 = 220.1 %TAS (ft/s)
q_bar = 49.6 %Dynamic Pressure (lbs/ft^2)
c_frac = 26.4 %CG Location
a_1 = 0 %Angle of Attack (Deg)

%Approach Flight Condition
h = 0
M = 0.096
U_1 = 107.1
q_bar = 13.6
c_frac = 26.4
a_1 = 4

%%%%%%%%%%%%%%%%Mass Data For All Flight Conditions%%%%%%%%%%%%%%%%%%%%%%
%Mass Data in Climb, Cruise, and Approach Conditions
W = 2650 %Weight (lbs)
m = W/g %Mass
I_xx_B = 948 %Moment of Inertia along the X-axis (slug-ft^2)
I_yy_B = 1346 %Moment of Inertia along the Y-axis (slug-ft^2)
I_zz_B = 1967 %Moment of Inertia along the Z-axis (slug-ft^2)
I_xz_B = 0 %Moment of Inertia along the XZ plane (slug-ft^2)

%%%%%%%%%%%%%%%%%%%%%%%Longitudinal Control and Hinge Momoment $$$$$$$$$
%Climb Longitudinal Control and Hinge Moment Derivatives
Ch_a = -0.0545
Ch_del_e = -0.594

%Cruise Longitudinal Control and Hinge Moment Derivatives
Ch_a = -0.0584
Ch_del_e = -0.585

%Approach Longitudinal Control and Hinge Moment Derivatives
Ch_a = -0.0549
Ch_del_e = -0.594

%%%%%%%%%%%%%%%%%%%Lateral-Directional Stability%%%%%%%%%%%%%%%%%%%%
%Climb Lateral-Directional Stabillity Derivatives
Cl_B = -0.0895
Cl_p = -0.487
Cl_r = 0.1869
Cy_B = -0.404
Cy_p = -0.145
Cy_r = 0.267
Cn_B = 0.0907
Cn_TB = 0
Cn_p = -0.0649
Cn_r = -0.1199

%Cruise Lateral-Directional Stability Derivatives
Cl_B = -0.0923
Cl_p = -0.484
Cl_r = 0.0798
Cy_B = -0.393
Cy_p = -0.075
Cy_r = 0.214
Cn_B = 0.0587
Cn_TB = 0
Cn_p = -0.0278
Cn_r = -0.0937

%Approach Lateral-Directional Stability Derivatives
Cl_B = -0.0969
Cl_p = -0.494
Cl_r = 0.2039
Cy_B = -0.303
Cy_p = -0.213
Cy_r = 0.201
Cn_B = 0.0701
Cn_TB = 0
Cn_p = -0.0960
Cn_r = -0.1151

%%%%%%%%%%%%%%%%%%%%Lateral-Directional Control and Hing Moment%%%%%%%%%
%Climb Lateral-Directional Control and Hinge Moment Derivatives
Cl_dela = 0.229
Cl_delr = 0.0147
Cy_dela = 0
Cy_delr = 0.187
Cn_dela = -0.0504
Cn_delr = -0.0805
Ch_aa = 0
Ch_dela = -0.369
Ch_Br = 0.0819
Ch_delr = -0.579

%Cruise Lateral-Directional Control and Hinge Moment Derivatives
Cl_dela = 0.229
Cl_delr = 0.0147
Cy_dela = 0
Cy_delr = 0.187
Cn_dela = -0.0216
Cn_delr = -0.0645
Ch_aa = 0
Ch_dela = -0.363
Ch_Br = 0.0819
Ch_delr = -0.567

%Approach Lateral-Directional Control and Hinge moment Derivatives
Cl_dela = 0.229
Cl_delr = 0.0147
Cy_dela = 0
Cy_delr = 0.187
Cn_dela = -0.0504
Cn_delr = -0.0805
Ch_aa = 0
Ch_dela = -0.369
Ch_Br = 0.0819
Ch_delr = -0.579

%Lateral-Directional Dimensional Stability Derivatives
YB = (q_bar*S*Cy_B)/m
Yp = (S*b*Cy_p)/(2*m*U_1)
Yr = (q_bar*S*b*Cy_r)/(2*m*U_1)
Y_dela = (q_bar*S*Cy_dela)/m
Y_delr = (q_bar*S*Cy_delr)/m
LB = (q_bar*S*b*Cl_B)/I_xx_B
Lp = (q_bar*S*(b^2)*Cl_p)/(2*I_xx_B*U_1)
Lr = (q_bar*S*(b^2)*Cl_r)/(2*I_xx_B*U_1)
L_dela = (q_bar*S*b*Cl_dela)/I_xx_B
L_delr = (q_bar*S*b*Cl_delr)/I_xx_B
NB = (q_bar*S*b*Cn_B)/I_zz_B
NT_B = (q_bar*S*b*Cn_TB)/I_zz_B
Np = (q_bar*S*(b^2)*Cn_p)/(2*I_zz_B*U_1)
Nr = (q_bar*S*(b^2)*Cn_r)/(2*I_zz_B*U_1)
N_dela = (q_bar*S*b*Cn_dela)/I_zz_B
N_delr = (q_bar*S*b*Cn_delr)/I_zz_B

%Perturbed Lateral-Directional Equations of Motion
A_1 = I_xz/I_xx_B
B_1 = I_xz/I_zz_B

%Denominator Polynomial
A2 = U_1(1-A_1*B_1)
B2 = -YB*(1-A_1*B_1)-U_1*(Lp+Nr+A_1*Np+B_1*Lr)
C2 = U_1*(Lp*Nr-Lr*Np)+YB*(Nr+Lp+A_1*Np+B_1*Lr)-Yp*(LB+NB*A_1+NT_B*A_1)...
    +U_1*(LB*B_1+NB+NT_B)-Yr*(LB*B_1+NB+NT_B)
D2 = -YB*(Lp*Nr-Lr*Np)+Yp*(LB*Nr-NB*Lr-NT_B*Lr)-g*cos(theta1)...
    *(LB+NB*A_1+NT_B*A_1)+U_1*(LB*Np-NB*Lp-NT_B*Lp)-Yr*(LB*Np-NB*Lp-NT_B*Lp)
E2 = g*cos(theta1)*(LB*Nr-NB*Lr-NT_B*Lr)

D_2 = [A2 B2 C2 D2 E2]

%Numerator Polynomial
AB = Y_dela*(1-A_1*B_1)+Y_delr*(1-A_1*B_1)
BB = -Y_dela*(Nr+Lp+A_1*Np+B_1*Lr)-Y_delr*(Nr+Lp+A_1*Np+B_1*Lr)+Yp*(L_dela+N_dela*A_1)+Yp*(L_delr+N_delr*A_1)...
    +Yr*(L_dela*B_1+N_dela)+Yr*(L_delr*B_1+N_delr)-U_1*(L_dela*B_1+N_dela)...
    +U_1*(L_delr*B_1+N_delr)
CB = Y_dela*(Lp*Nr-Np*Lr)+Y_delr*(Lp*Nr-Np*Lr)+Yp*(N_dela*Lr-L_dela*Nr)+Yp*(N_delr*Lr-L_delr*Nr)...
    +g*cos(theta1)*(L_dela+N_dela*A_1)+g*cos(theta1)*(L_delr+N_delr*A_1)...
    +Yr*(L_dela*Np-N_dela*Lp)+Yr*(L_delr*Np-N_delr*Lp)+U_1*(L_dela*Np-N_dela*Lp)+U_1*(L_delr*Np-N_delr*Lp)
DB = g*cos(theta1)*(N_dela*Lr-L_dela*Nr)+g*cos(theta1)*(N_delr*Lr-L_delr*Nr)
N_B = [AB BB CB DB]
%Transfer Function
sys = tf(N_B,D_2)
damp(sys)

Answer 1

[Sam Chak]

Hi @Jeffrey Lewis

There are a few ways to do it. But since you have written the data into the code, I think probably the simplest way is to write a function file for eash flight condition, and then call the specific function. The example is shown below:

1. Calling climbsys

climbsys

Climb_sys =
 
     8.375 s^3 - 200.8 s^2 - 1.251e04 s + 195.1
  ------------------------------------------------
  133.5 s^4 + 1392 s^3 + 2612 s^2 + 9028 s - 255.2
 
Continuous-time transfer function.

                                                                       
         Pole              Damping       Frequency      Time Constant  
                                       (rad/seconds)      (seconds)    
                                                                       
  2.80e-02                -1.00e+00       2.80e-02        -3.57e+01    
 -6.80e-01 + 2.65e+00i     2.48e-01       2.74e+00         1.47e+00    
 -6.80e-01 - 2.65e+00i     2.48e-01       2.74e+00         1.47e+00    
 -9.09e+00                 1.00e+00       9.09e+00         1.10e-01    

2. Calling cruisesys

cruisesys

Cruise_sys =
 
       19.59 s^3 - 1240 s^2 - 4.263e04 s + 2174
  --------------------------------------------------
  220.1 s^4 + 3163 s^3 + 6232 s^2 + 3.028e04 s + 540
 
Continuous-time transfer function.

                                                                       
         Pole              Damping       Frequency      Time Constant  
                                       (rad/seconds)      (seconds)    
                                                                       
 -1.79e-02                 1.00e+00       1.79e-02         5.59e+01    
 -6.74e-01 + 3.18e+00i     2.08e-01       3.25e+00         1.48e+00    
 -6.74e-01 - 3.18e+00i     2.08e-01       3.25e+00         1.48e+00    
 -1.30e+01                 1.00e+00       1.30e+01         7.69e-02    

3. Calling approachsys

approachsys

Approach_sys =
 
       5.373 s^3 - 102.2 s^2 - 5475 s + 28.61
  -------------------------------------------------
  107.1 s^4 + 897.5 s^3 + 1294 s^2 + 3403 s - 66.09
 
Continuous-time transfer function.

                                                                       
         Pole              Damping       Frequency      Time Constant  
                                       (rad/seconds)      (seconds)    
                                                                       
  1.93e-02                -1.00e+00       1.93e-02        -5.19e+01    
 -5.37e-01 + 2.02e+00i     2.57e-01       2.09e+00         1.86e+00    
 -5.37e-01 - 2.02e+00i     2.57e-01       2.09e+00         1.86e+00    
 -7.32e+00                 1.00e+00       7.32e+00         1.37e-01    

function climbsys

    g      = 32.174;
    theta1 = 0;
    
    % Cessna design parameters
    S = 174;    % Planform Area (ft^2)
    c = 4.9;    % Chord length (ft)
    b = 36.0;   % Wing span (ft)
    
    %%%%%%%%%%%%%%%%Flight Conditions%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Climb Flight Condition
    altitude = 0;
    M        = 0.120;
    U_1      = 133.5;
    q_bar    = 21.2;
    c_frac   = 26.2;
    a_1      = 5.4;
    
    % % Cruise Flight Condition
    % h      = 5000;      % Altitude (ft)
    % M      = 0.201;     % Mach Number
    % U_1    = 220.1;     % TAS (ft/s)
    % q_bar  = 49.6;      % Dynamic Pressure (lbs/ft^2)
    % c_frac = 26.4;      % CG Location
    % a_1    = 0;         % Angle of Attack (Deg)
    % 
    % % Approach Flight Condition
    % h      = 0;
    % M      = 0.096;
    % U_1    = 107.1;
    % q_bar  = 13.6;
    % c_frac = 26.4;
    % a_1    = 4;
    
    %%%%%%%%%%%%%%%%Mass Data For All Flight Conditions%%%%%%%%%%%%%%%%%%%%%%
    % Mass Data in Climb, Cruise, and Approach Conditions
    W      = 2650;      % Weight (lbs)
    m      = W/g;       % Mass
    I_xx_B = 948;       % Moment of Inertia along the X-axis (slug-ft^2)
    I_yy_B = 1346;      % Moment of Inertia along the Y-axis (slug-ft^2)
    I_zz_B = 1967;      % Moment of Inertia along the Z-axis (slug-ft^2)
    I_xz_B = 0;         % Moment of Inertia along the XZ plane (slug-ft^2)
    
    %%%%%%%%%%%%%%%%%%%%%%%Longitudinal Control and Hinge Momoment $$$$$$$$$
    % Climb Longitudinal Control and Hinge Moment Derivatives
    Ch_a     = -0.0545;
    Ch_del_e = -0.594;
    
    % % Cruise Longitudinal Control and Hinge Moment Derivatives
    % Ch_a     = -0.0584;
    % Ch_del_e = -0.585;
    % 
    % % Approach Longitudinal Control and Hinge Moment Derivatives
    % Ch_a     = -0.0549;
    % Ch_del_e = -0.594;
    
    %%%%%%%%%%%%%%%%%%%Lateral-Directional Stability%%%%%%%%%%%%%%%%%%%%
    % Climb Lateral-Directional Stabillity Derivatives
    Cl_B  = -0.0895;
    Cl_p  = -0.487;
    Cl_r  =  0.1869;
    Cy_B  = -0.404;
    Cy_p  = -0.145;
    Cy_r  =  0.267;
    Cn_B  =  0.0907;
    Cn_TB =  0;
    Cn_p  = -0.0649;
    Cn_r  = -0.1199;
    
    % % Cruise Lateral-Directional Stability Derivatives
    % Cl_B  = -0.0923;
    % Cl_p  = -0.484;
    % Cl_r  =  0.0798;
    % Cy_B  = -0.393;
    % Cy_p  = -0.075;
    % Cy_r  =  0.214;
    % Cn_B  =  0.0587;
    % Cn_TB =  0;
    % Cn_p  = -0.0278;
    % Cn_r  = -0.0937;
    % 
    % % Approach Lateral-Directional Stability Derivatives
    % Cl_B  = -0.0969;
    % Cl_p  = -0.494;
    % Cl_r  =  0.2039;
    % Cy_B  = -0.303;
    % Cy_p  = -0.213;
    % Cy_r  =  0.201;
    % Cn_B  =  0.0701;
    % Cn_TB =  0;
    % Cn_p  = -0.0960;
    % Cn_r  = -0.1151;
    
    %%%%%%%%%%%%%%%%%%%%Lateral-Directional Control and Hing Moment%%%%%%%%%
    % Climb Lateral-Directional Control and Hinge Moment Derivatives
    Cl_dela =  0.229;
    Cl_delr =  0.0147;
    Cy_dela =  0;
    Cy_delr =  0.187;
    Cn_dela = -0.0504;
    Cn_delr = -0.0805;
    Ch_aa   =  0;
    Ch_dela = -0.369;
    Ch_Br   =  0.0819;
    Ch_delr = -0.579;
    
    % % Cruise Lateral-Directional Control and Hinge Moment Derivatives
    % Cl_dela =  0.229;
    % Cl_delr =  0.0147;
    % Cy_dela =  0;
    % Cy_delr =  0.187;
    % Cn_dela = -0.0216;
    % Cn_delr = -0.0645;
    % Ch_aa   =  0;
    % Ch_dela = -0.363;
    % Ch_Br   =  0.0819;
    % Ch_delr = -0.567;
    % 
    % % Approach Lateral-Directional Control and Hinge moment Derivatives
    % Cl_dela =  0.229;
    % Cl_delr =  0.0147;
    % Cy_dela =  0;
    % Cy_delr =  0.187;
    % Cn_dela = -0.0504;
    % Cn_delr = -0.0805;
    % Ch_aa   =  0;
    % Ch_dela = -0.369;
    % Ch_Br   =  0.0819;
    % Ch_delr = -0.579;
    
    % Lateral-Directional Dimensional Stability Derivatives
    YB     = (q_bar*S*Cy_B)/m;
    Yp     = (S*b*Cy_p)/(2*m*U_1);
    Yr     = (q_bar*S*b*Cy_r)/(2*m*U_1);
    Y_dela = (q_bar*S*Cy_dela)/m;
    Y_delr = (q_bar*S*Cy_delr)/m;
    LB     = (q_bar*S*b*Cl_B)/I_xx_B;
    Lp     = (q_bar*S*(b^2)*Cl_p)/(2*I_xx_B*U_1);
    Lr     = (q_bar*S*(b^2)*Cl_r)/(2*I_xx_B*U_1);
    L_dela = (q_bar*S*b*Cl_dela)/I_xx_B;
    L_delr = (q_bar*S*b*Cl_delr)/I_xx_B;
    NB     = (q_bar*S*b*Cn_B)/I_zz_B;
    NT_B   = (q_bar*S*b*Cn_TB)/I_zz_B;
    Np     = (q_bar*S*(b^2)*Cn_p)/(2*I_zz_B*U_1);
    Nr     = (q_bar*S*(b^2)*Cn_r)/(2*I_zz_B*U_1);
    N_dela = (q_bar*S*b*Cn_dela)/I_zz_B;
    N_delr = (q_bar*S*b*Cn_delr)/I_zz_B;
    
    % Perturbed Lateral-Directional Equations of Motion
    A_1 = I_xz_B/I_xx_B;
    B_1 = I_xz_B/I_zz_B;
    
    % Numerator Polynomial
    AB  =   Y_dela*(1 - A_1*B_1) + Y_delr*(1 - A_1*B_1);
    BB  = - Y_dela*(Nr + Lp + A_1*Np + B_1*Lr) - Y_delr*(Nr + Lp + A_1*Np + B_1*Lr) + Yp*(L_dela + N_dela*A_1) + Yp*(L_delr + N_delr*A_1) + Yr*(L_dela*B_1 + N_dela) + Yr*(L_delr*B_1 + N_delr) - U_1*(L_dela*B_1 + N_dela) + U_1*(L_delr*B_1 + N_delr);
    CB  =   Y_dela*(Lp*Nr - Np*Lr) + Y_delr*(Lp*Nr - Np*Lr) + Yp*(N_dela*Lr - L_dela*Nr) + Yp*(N_delr*Lr - L_delr*Nr) + g*cos(theta1)*(L_dela + N_dela*A_1) + g*cos(theta1)*(L_delr + N_delr*A_1) + Yr*(L_dela*Np - N_dela*Lp) + Yr*(L_delr*Np - N_delr*Lp) + U_1*(L_dela*Np - N_dela*Lp) + U_1*(L_delr*Np - N_delr*Lp);
    DB  =   g*cos(theta1)*(N_dela*Lr-L_dela*Nr)+g*cos(theta1)*(N_delr*Lr-L_delr*Nr);
    N_B =   [AB BB CB DB];
    
    % Denominator Polynomial
    A2  =   U_1*(1 - A_1*B_1);
    B2  = - YB*(1 - A_1*B_1) - U_1*(Lp + Nr + A_1*Np + B_1*Lr);
    C2  =   U_1*(Lp*Nr - Lr*Np) + YB*(Nr + Lp + A_1*Np + B_1*Lr) - Yp*(LB + NB*A_1 + NT_B*A_1) + U_1*(LB*B_1 + NB + NT_B) - Yr*(LB*B_1 + NB + NT_B);
    D2  = - YB*(Lp*Nr - Lr*Np) + Yp*(LB*Nr - NB*Lr - NT_B*Lr) - g*cos(theta1)*(LB + NB*A_1 + NT_B*A_1) + U_1*(LB*Np - NB*Lp - NT_B*Lp) - Yr*(LB*Np - NB*Lp - NT_B*Lp);
    E2  =   g*cos(theta1)*(LB*Nr - NB*Lr - NT_B*Lr);
    D_2 =   [A2 B2 C2 D2 E2];
    
    % Transfer Function
    Climb_sys = tf(N_B, D_2)
    damp(Climb_sys)
    step(Climb_sys)

end

function cruisesys

    g      = 32.174;
    theta1 = 0;
    
    % Cessna design parameters
    S = 174;    % Planform Area (ft^2)
    c = 4.9;    % Chord length (ft)
    b = 36.0;   % Wing span (ft)
    
    %%%%%%%%%%%%%%%%Flight Conditions%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % % Climb Flight Condition
    % altitude = 0;
    % M        = 0.120;
    % U_1      = 133.5;
    % q_bar    = 21.2;
    % c_frac   = 26.2;
    % a_1      = 5.4;
    
    % Cruise Flight Condition
    h      = 5000;      % Altitude (ft)
    M      = 0.201;     % Mach Number
    U_1    = 220.1;     % TAS (ft/s)
    q_bar  = 49.6;      % Dynamic Pressure (lbs/ft^2)
    c_frac = 26.4;      % CG Location
    a_1    = 0;         % Angle of Attack (Deg)
    
    % % Approach Flight Condition
    % h      = 0;
    % M      = 0.096;
    % U_1    = 107.1;
    % q_bar  = 13.6;
    % c_frac = 26.4;
    % a_1    = 4;
    
    %%%%%%%%%%%%%%%%Mass Data For All Flight Conditions%%%%%%%%%%%%%%%%%%%%%%
    % Mass Data in Climb, Cruise, and Approach Conditions
    W      = 2650;      % Weight (lbs)
    m      = W/g;       % Mass
    I_xx_B = 948;       % Moment of Inertia along the X-axis (slug-ft^2)
    I_yy_B = 1346;      % Moment of Inertia along the Y-axis (slug-ft^2)
    I_zz_B = 1967;      % Moment of Inertia along the Z-axis (slug-ft^2)
    I_xz_B = 0;         % Moment of Inertia along the XZ plane (slug-ft^2)
    
    %%%%%%%%%%%%%%%%%%%%%%%Longitudinal Control and Hinge Momoment $$$$$$$$$
    % % Climb Longitudinal Control and Hinge Moment Derivatives
    % Ch_a     = -0.0545;
    % Ch_del_e = -0.594;
    
    % Cruise Longitudinal Control and Hinge Moment Derivatives
    Ch_a     = -0.0584;
    Ch_del_e = -0.585;
    
    % % Approach Longitudinal Control and Hinge Moment Derivatives
    % Ch_a     = -0.0549;
    % Ch_del_e = -0.594;
    
    %%%%%%%%%%%%%%%%%%%Lateral-Directional Stability%%%%%%%%%%%%%%%%%%%%
    % % Climb Lateral-Directional Stabillity Derivatives
    % Cl_B  = -0.0895;
    % Cl_p  = -0.487;
    % Cl_r  =  0.1869;
    % Cy_B  = -0.404;
    % Cy_p  = -0.145;
    % Cy_r  =  0.267;
    % Cn_B  =  0.0907;
    % Cn_TB =  0;
    % Cn_p  = -0.0649;
    % Cn_r  = -0.1199;
    
    % Cruise Lateral-Directional Stability Derivatives
    Cl_B  = -0.0923;
    Cl_p  = -0.484;
    Cl_r  =  0.0798;
    Cy_B  = -0.393;
    Cy_p  = -0.075;
    Cy_r  =  0.214;
    Cn_B  =  0.0587;
    Cn_TB =  0;
    Cn_p  = -0.0278;
    Cn_r  = -0.0937;
    
    % % Approach Lateral-Directional Stability Derivatives
    % Cl_B  = -0.0969;
    % Cl_p  = -0.494;
    % Cl_r  =  0.2039;
    % Cy_B  = -0.303;
    % Cy_p  = -0.213;
    % Cy_r  =  0.201;
    % Cn_B  =  0.0701;
    % Cn_TB =  0;
    % Cn_p  = -0.0960;
    % Cn_r  = -0.1151;
    
    %%%%%%%%%%%%%%%%%%%%Lateral-Directional Control and Hing Moment%%%%%%%%%
    % % Climb Lateral-Directional Control and Hinge Moment Derivatives
    % Cl_dela =  0.229;
    % Cl_delr =  0.0147;
    % Cy_dela =  0;
    % Cy_delr =  0.187;
    % Cn_dela = -0.0504;
    % Cn_delr = -0.0805;
    % Ch_aa   =  0;
    % Ch_dela = -0.369;
    % Ch_Br   =  0.0819;
    % Ch_delr = -0.579;
    
    % Cruise Lateral-Directional Control and Hinge Moment Derivatives
    Cl_dela =  0.229;
    Cl_delr =  0.0147;
    Cy_dela =  0;
    Cy_delr =  0.187;
    Cn_dela = -0.0216;
    Cn_delr = -0.0645;
    Ch_aa   =  0;
    Ch_dela = -0.363;
    Ch_Br   =  0.0819;
    Ch_delr = -0.567;
    
    % % Approach Lateral-Directional Control and Hinge moment Derivatives
    % Cl_dela =  0.229;
    % Cl_delr =  0.0147;
    % Cy_dela =  0;
    % Cy_delr =  0.187;
    % Cn_dela = -0.0504;
    % Cn_delr = -0.0805;
    % Ch_aa   =  0;
    % Ch_dela = -0.369;
    % Ch_Br   =  0.0819;
    % Ch_delr = -0.579;
    
    % Lateral-Directional Dimensional Stability Derivatives
    YB     = (q_bar*S*Cy_B)/m;
    Yp     = (S*b*Cy_p)/(2*m*U_1);
    Yr     = (q_bar*S*b*Cy_r)/(2*m*U_1);
    Y_dela = (q_bar*S*Cy_dela)/m;
    Y_delr = (q_bar*S*Cy_delr)/m;
    LB     = (q_bar*S*b*Cl_B)/I_xx_B;
    Lp     = (q_bar*S*(b^2)*Cl_p)/(2*I_xx_B*U_1);
    Lr     = (q_bar*S*(b^2)*Cl_r)/(2*I_xx_B*U_1);
    L_dela = (q_bar*S*b*Cl_dela)/I_xx_B;
    L_delr = (q_bar*S*b*Cl_delr)/I_xx_B;
    NB     = (q_bar*S*b*Cn_B)/I_zz_B;
    NT_B   = (q_bar*S*b*Cn_TB)/I_zz_B;
    Np     = (q_bar*S*(b^2)*Cn_p)/(2*I_zz_B*U_1);
    Nr     = (q_bar*S*(b^2)*Cn_r)/(2*I_zz_B*U_1);
    N_dela = (q_bar*S*b*Cn_dela)/I_zz_B;
    N_delr = (q_bar*S*b*Cn_delr)/I_zz_B;
    
    % Perturbed Lateral-Directional Equations of Motion
    A_1 = I_xz_B/I_xx_B;
    B_1 = I_xz_B/I_zz_B;
    
    % Numerator Polynomial
    AB  =   Y_dela*(1 - A_1*B_1) + Y_delr*(1 - A_1*B_1);
    BB  = - Y_dela*(Nr + Lp + A_1*Np + B_1*Lr) - Y_delr*(Nr + Lp + A_1*Np + B_1*Lr) + Yp*(L_dela + N_dela*A_1) + Yp*(L_delr + N_delr*A_1) + Yr*(L_dela*B_1 + N_dela) + Yr*(L_delr*B_1 + N_delr) - U_1*(L_dela*B_1 + N_dela) + U_1*(L_delr*B_1 + N_delr);
    CB  =   Y_dela*(Lp*Nr - Np*Lr) + Y_delr*(Lp*Nr - Np*Lr) + Yp*(N_dela*Lr - L_dela*Nr) + Yp*(N_delr*Lr - L_delr*Nr) + g*cos(theta1)*(L_dela + N_dela*A_1) + g*cos(theta1)*(L_delr + N_delr*A_1) + Yr*(L_dela*Np - N_dela*Lp) + Yr*(L_delr*Np - N_delr*Lp) + U_1*(L_dela*Np - N_dela*Lp) + U_1*(L_delr*Np - N_delr*Lp);
    DB  =   g*cos(theta1)*(N_dela*Lr-L_dela*Nr)+g*cos(theta1)*(N_delr*Lr-L_delr*Nr);
    N_B =   [AB BB CB DB];
    
    % Denominator Polynomial
    A2  =   U_1*(1 - A_1*B_1);
    B2  = - YB*(1 - A_1*B_1) - U_1*(Lp + Nr + A_1*Np + B_1*Lr);
    C2  =   U_1*(Lp*Nr - Lr*Np) + YB*(Nr + Lp + A_1*Np + B_1*Lr) - Yp*(LB + NB*A_1 + NT_B*A_1) + U_1*(LB*B_1 + NB + NT_B) - Yr*(LB*B_1 + NB + NT_B);
    D2  = - YB*(Lp*Nr - Lr*Np) + Yp*(LB*Nr - NB*Lr - NT_B*Lr) - g*cos(theta1)*(LB + NB*A_1 + NT_B*A_1) + U_1*(LB*Np - NB*Lp - NT_B*Lp) - Yr*(LB*Np - NB*Lp - NT_B*Lp);
    E2  =   g*cos(theta1)*(LB*Nr - NB*Lr - NT_B*Lr);
    D_2 =   [A2 B2 C2 D2 E2];
    
    % Transfer Function
    Cruise_sys = tf(N_B, D_2)
    damp(Cruise_sys)
    step(Cruise_sys)

end

function approachsys

    g      = 32.174;
    theta1 = 0;
    
    % Cessna design parameters
    S = 174;    % Planform Area (ft^2)
    c = 4.9;    % Chord length (ft)
    b = 36.0;   % Wing span (ft)
    
    %%%%%%%%%%%%%%%%Flight Conditions%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % % Climb Flight Condition
    % altitude = 0;
    % M        = 0.120;
    % U_1      = 133.5;
    % q_bar    = 21.2;
    % c_frac   = 26.2;
    % a_1      = 5.4;
    % 
    % % Cruise Flight Condition
    % h      = 5000;      % Altitude (ft)
    % M      = 0.201;     % Mach Number
    % U_1    = 220.1;     % TAS (ft/s)
    % q_bar  = 49.6;      % Dynamic Pressure (lbs/ft^2)
    % c_frac = 26.4;      % CG Location
    % a_1    = 0;         % Angle of Attack (Deg)
    
    % Approach Flight Condition
    h      = 0;
    M      = 0.096;
    U_1    = 107.1;
    q_bar  = 13.6;
    c_frac = 26.4;
    a_1    = 4;
    
    %%%%%%%%%%%%%%%%Mass Data For All Flight Conditions%%%%%%%%%%%%%%%%%%%%%%
    % Mass Data in Climb, Cruise, and Approach Conditions
    W      = 2650;      % Weight (lbs)
    m      = W/g;       % Mass
    I_xx_B = 948;       % Moment of Inertia along the X-axis (slug-ft^2)
    I_yy_B = 1346;      % Moment of Inertia along the Y-axis (slug-ft^2)
    I_zz_B = 1967;      % Moment of Inertia along the Z-axis (slug-ft^2)
    I_xz_B = 0;         % Moment of Inertia along the XZ plane (slug-ft^2)
    
    %%%%%%%%%%%%%%%%%%%%%%%Longitudinal Control and Hinge Momoment $$$$$$$$$
    % % Climb Longitudinal Control and Hinge Moment Derivatives
    % Ch_a     = -0.0545;
    % Ch_del_e = -0.594;
    % 
    % % Cruise Longitudinal Control and Hinge Moment Derivatives
    % Ch_a     = -0.0584;
    % Ch_del_e = -0.585;
    
    % Approach Longitudinal Control and Hinge Moment Derivatives
    Ch_a     = -0.0549;
    Ch_del_e = -0.594;
    
    %%%%%%%%%%%%%%%%%%%Lateral-Directional Stability%%%%%%%%%%%%%%%%%%%%
    % % Climb Lateral-Directional Stabillity Derivatives
    % Cl_B  = -0.0895;
    % Cl_p  = -0.487;
    % Cl_r  =  0.1869;
    % Cy_B  = -0.404;
    % Cy_p  = -0.145;
    % Cy_r  =  0.267;
    % Cn_B  =  0.0907;
    % Cn_TB =  0;
    % Cn_p  = -0.0649;
    % Cn_r  = -0.1199;
    % 
    % % Cruise Lateral-Directional Stability Derivatives
    % Cl_B  = -0.0923;
    % Cl_p  = -0.484;
    % Cl_r  =  0.0798;
    % Cy_B  = -0.393;
    % Cy_p  = -0.075;
    % Cy_r  =  0.214;
    % Cn_B  =  0.0587;
    % Cn_TB =  0;
    % Cn_p  = -0.0278;
    % Cn_r  = -0.0937;
    
    % Approach Lateral-Directional Stability Derivatives
    Cl_B  = -0.0969;
    Cl_p  = -0.494;
    Cl_r  =  0.2039;
    Cy_B  = -0.303;
    Cy_p  = -0.213;
    Cy_r  =  0.201;
    Cn_B  =  0.0701;
    Cn_TB =  0;
    Cn_p  = -0.0960;
    Cn_r  = -0.1151;
    
    %%%%%%%%%%%%%%%%%%%%Lateral-Directional Control and Hing Moment%%%%%%%%%
    % % Climb Lateral-Directional Control and Hinge Moment Derivatives
    % Cl_dela =  0.229;
    % Cl_delr =  0.0147;
    % Cy_dela =  0;
    % Cy_delr =  0.187;
    % Cn_dela = -0.0504;
    % Cn_delr = -0.0805;
    % Ch_aa   =  0;
    % Ch_dela = -0.369;
    % Ch_Br   =  0.0819;
    % Ch_delr = -0.579;
    % 
    % % Cruise Lateral-Directional Control and Hinge Moment Derivatives
    % Cl_dela =  0.229;
    % Cl_delr =  0.0147;
    % Cy_dela =  0;
    % Cy_delr =  0.187;
    % Cn_dela = -0.0216;
    % Cn_delr = -0.0645;
    % Ch_aa   =  0;
    % Ch_dela = -0.363;
    % Ch_Br   =  0.0819;
    % Ch_delr = -0.567;
    
    % Approach Lateral-Directional Control and Hinge moment Derivatives
    Cl_dela =  0.229;
    Cl_delr =  0.0147;
    Cy_dela =  0;
    Cy_delr =  0.187;
    Cn_dela = -0.0504;
    Cn_delr = -0.0805;
    Ch_aa   =  0;
    Ch_dela = -0.369;
    Ch_Br   =  0.0819;
    Ch_delr = -0.579;
    
    % Lateral-Directional Dimensional Stability Derivatives
    YB     = (q_bar*S*Cy_B)/m;
    Yp     = (S*b*Cy_p)/(2*m*U_1);
    Yr     = (q_bar*S*b*Cy_r)/(2*m*U_1);
    Y_dela = (q_bar*S*Cy_dela)/m;
    Y_delr = (q_bar*S*Cy_delr)/m;
    LB     = (q_bar*S*b*Cl_B)/I_xx_B;
    Lp     = (q_bar*S*(b^2)*Cl_p)/(2*I_xx_B*U_1);
    Lr     = (q_bar*S*(b^2)*Cl_r)/(2*I_xx_B*U_1);
    L_dela = (q_bar*S*b*Cl_dela)/I_xx_B;
    L_delr = (q_bar*S*b*Cl_delr)/I_xx_B;
    NB     = (q_bar*S*b*Cn_B)/I_zz_B;
    NT_B   = (q_bar*S*b*Cn_TB)/I_zz_B;
    Np     = (q_bar*S*(b^2)*Cn_p)/(2*I_zz_B*U_1);
    Nr     = (q_bar*S*(b^2)*Cn_r)/(2*I_zz_B*U_1);
    N_dela = (q_bar*S*b*Cn_dela)/I_zz_B;
    N_delr = (q_bar*S*b*Cn_delr)/I_zz_B;
    
    % Perturbed Lateral-Directional Equations of Motion
    A_1 = I_xz_B/I_xx_B;
    B_1 = I_xz_B/I_zz_B;
    
    % Numerator Polynomial
    AB  =   Y_dela*(1 - A_1*B_1) + Y_delr*(1 - A_1*B_1);
    BB  = - Y_dela*(Nr + Lp + A_1*Np + B_1*Lr) - Y_delr*(Nr + Lp + A_1*Np + B_1*Lr) + Yp*(L_dela + N_dela*A_1) + Yp*(L_delr + N_delr*A_1) + Yr*(L_dela*B_1 + N_dela) + Yr*(L_delr*B_1 + N_delr) - U_1*(L_dela*B_1 + N_dela) + U_1*(L_delr*B_1 + N_delr);
    CB  =   Y_dela*(Lp*Nr - Np*Lr) + Y_delr*(Lp*Nr - Np*Lr) + Yp*(N_dela*Lr - L_dela*Nr) + Yp*(N_delr*Lr - L_delr*Nr) + g*cos(theta1)*(L_dela + N_dela*A_1) + g*cos(theta1)*(L_delr + N_delr*A_1) + Yr*(L_dela*Np - N_dela*Lp) + Yr*(L_delr*Np - N_delr*Lp) + U_1*(L_dela*Np - N_dela*Lp) + U_1*(L_delr*Np - N_delr*Lp);
    DB  =   g*cos(theta1)*(N_dela*Lr-L_dela*Nr)+g*cos(theta1)*(N_delr*Lr-L_delr*Nr);
    N_B =   [AB BB CB DB];
    
    % Denominator Polynomial
    A2  =   U_1*(1 - A_1*B_1);
    B2  = - YB*(1 - A_1*B_1) - U_1*(Lp + Nr + A_1*Np + B_1*Lr);
    C2  =   U_1*(Lp*Nr - Lr*Np) + YB*(Nr + Lp + A_1*Np + B_1*Lr) - Yp*(LB + NB*A_1 + NT_B*A_1) + U_1*(LB*B_1 + NB + NT_B) - Yr*(LB*B_1 + NB + NT_B);
    D2  = - YB*(Lp*Nr - Lr*Np) + Yp*(LB*Nr - NB*Lr - NT_B*Lr) - g*cos(theta1)*(LB + NB*A_1 + NT_B*A_1) + U_1*(LB*Np - NB*Lp - NT_B*Lp) - Yr*(LB*Np - NB*Lp - NT_B*Lp);
    E2  =   g*cos(theta1)*(LB*Nr - NB*Lr - NT_B*Lr);
    D_2 =   [A2 B2 C2 D2 E2];
    
    % Transfer Function
    Approach_sys = tf(N_B, D_2)
    damp(Approach_sys)
    step(Approach_sys)

end

[Jeffrey Lewis]

This is the approach I was looking for, thank you! The only issue I am facing with this structure is I am not getting an output on my end. Is there anything extra I need to do at the end to call the functions or should I expect the out put to show my TF as you have displayed?

I have written code for a previous homework that only includes one flight condition and it still outputs the TF as you have shown but now that I have all three conditions in this function structure, MATLAB is not outputting anything.

clc 
clear all

function climb_sys

    g = 32.174;
    theta1 = 0;
    %Cessna design parameters
    S = 174; %Planform Area (ft^2)
    c = 4.9; %Chord length (ft)
    b = 36.0; %Wing span (ft)

    %%%%%%%%%%%%%%%%Flight Conditions%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Climb Flight Condition
    altitude = 0;
    M = 0.120;
    U_1 = 133.5;
    q_bar = 21.2;
    c_frac = 26.2;
    a_1 = 5.4;

    %%%%%%%%%%%%%%%%Mass Data For All Flight Conditions%%%%%%%%%%%%%%%%%%%%%%
    %Mass Data in Climb, Cruise, and Approach Conditions
    W = 2650; %Weight (lbs)
    m = W/g; %Mass
    I_xx_B = 948; %Moment of Inertia along the X-axis (slug-ft^2)
    I_yy_B = 1346; %Moment of Inertia along the Y-axis (slug-ft^2)
    I_zz_B = 1967; %Moment of Inertia along the Z-axis (slug-ft^2)
    I_xz_B = 0; %Moment of Inertia along the XZ plane (slug-ft^2)

    %%%%%%%%%%%%%%%%%%%%%%%Longitudinal Control and Hinge Momoment $$$$$$$$$
    %Climb Longitudinal Control and Hinge Moment Derivatives
    Ch_a = -0.0545;
    Ch_del_e = -0.594;

    %%%%%%%%%%%%%%%%%%%Lateral-Directional Stability%%%%%%%%%%%%%%%%%%%%
    %Climb Lateral-Directional Stabillity Derivatives
    Cl_B = -0.0895;
    Cl_p = -0.487;
    Cl_r = 0.1869;
    Cy_B = -0.404;
    Cy_p = -0.145;
    Cy_r = 0.267;
    Cn_B = 0.0907;
    Cn_TB = 0;
    Cn_p = -0.0649;
    Cn_r = -0.1199;

    %%%%%%%%%%%%%%%%%%%%Lateral-Directional Control and Hing Moment%%%%%%%%%
    %Climb Lateral-Directional Control and Hinge Moment Derivatives
    Cl_dela = 0.229;
    Cl_delr = 0.0147;
    Cy_dela = 0;
    Cy_delr = 0.187;
    Cn_dela = -0.0504;
    Cn_delr = -0.0805;
    Ch_aa = 0;
    Ch_dela = -0.369;
    Ch_Br = 0.0819;
    Ch_delr = -0.579;

    %Lateral-Directional Dimensional Stability Derivatives
    YB = (q_bar*S*Cy_B)/m;
    Yp = (S*b*Cy_p)/(2*m*U_1);
    Yr = (q_bar*S*b*Cy_r)/(2*m*U_1);
    Y_dela = (q_bar*S*Cy_dela)/m;
    Y_delr = (q_bar*S*Cy_delr)/m;
    LB = (q_bar*S*b*Cl_B)/I_xx_B;
    Lp = (q_bar*S*(b^2)*Cl_p)/(2*I_xx_B*U_1);
    Lr = (q_bar*S*(b^2)*Cl_r)/(2*I_xx_B*U_1);
    L_dela = (q_bar*S*b*Cl_dela)/I_xx_B;
    L_delr = (q_bar*S*b*Cl_delr)/I_xx_B;
    NB = (q_bar*S*b*Cn_B)/I_zz_B;
    NT_B = (q_bar*S*b*Cn_TB)/I_zz_B;
    Np = (q_bar*S*(b^2)*Cn_p)/(2*I_zz_B*U_1);
    Nr = (q_bar*S*(b^2)*Cn_r)/(2*I_zz_B*U_1);
    N_dela = (q_bar*S*b*Cn_dela)/I_zz_B;
    N_delr = (q_bar*S*b*Cn_delr)/I_zz_B;

    %Perturbed Lateral-Directional Equations of Motion
    A_1 = I_xz/I_xx_B;
    B_1 = I_xz/I_zz_B;

    %Denominator Polynomial
    A2 = U_1(1-A_1*B_1);
    B2 = -YB*(1-A_1*B_1)-U_1*(Lp+Nr+A_1*Np+B_1*Lr);
    C2 = U_1*(Lp*Nr-Lr*Np)+YB*(Nr+Lp+A_1*Np+B_1*Lr)-Yp*(LB+NB*A_1+NT_B*A_1)...
    +U_1*(LB*B_1+NB+NT_B)-Yr*(LB*B_1+NB+NT_B);
    D2 = -YB*(Lp*Nr-Lr*Np)+Yp*(LB*Nr-NB*Lr-NT_B*Lr)-g*cos(theta1)...
    *(LB+NB*A_1+NT_B*A_1)+U_1*(LB*Np-NB*Lp-NT_B*Lp)-Yr*(LB*Np-NB*Lp-NT_B*Lp);
    E2 = g*cos(theta1)*(LB*Nr-NB*Lr-NT_B*Lr);

    D_2 = [A2 B2 C2 D2 E2];

    %Numerator Polynomial
    AB = Y_dela*(1-A_1*B_1)+Y_delr*(1-A_1*B_1);
    BB = -Y_dela*(Nr+Lp+A_1*Np+B_1*Lr)-Y_delr*(Nr+Lp+A_1*Np+B_1*Lr)+Yp*(L_dela+N_dela*A_1)+Yp*(L_delr+N_delr*A_1)...
        +Yr*(L_dela*B_1+N_dela)+Yr*(L_delr*B_1+N_delr)-U_1*(L_dela*B_1+N_dela)...
        +U_1*(L_delr*B_1+N_delr);
    CB = Y_dela*(Lp*Nr-Np*Lr)+Y_delr*(Lp*Nr-Np*Lr)+Yp*(N_dela*Lr-L_dela*Nr)+Yp*(N_delr*Lr-L_delr*Nr)...
        +g*cos(theta1)*(L_dela+N_dela*A_1)+g*cos(theta1)*(L_delr+N_delr*A_1)...
        +Yr*(L_dela*Np-N_dela*Lp)+Yr*(L_delr*Np-N_delr*Lp)+U_1*(L_dela*Np-N_dela*Lp)+U_1*(L_delr*Np-N_delr*Lp);
    DB = g*cos(theta1)*(N_dela*Lr-L_dela*Nr)+g*cos(theta1)*(N_delr*Lr-L_delr*Nr);
    N_B = [AB BB CB DB];

    %Transfer Function
    climb_sys = tf(N_B,D_2);
    damp(climb_sys);
end
function cruise_sys
    g = 32.174
    theta1 = 0
    %Cessna design parameters
    S = 174 %Planform Area (ft^2)
    c = 4.9 %Chord length (ft)
    b = 36.0 %Wing span (ft)

    %%%%%%%%%%%%%%%%Flight Conditions%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %Cruise Flight Condition
    h = 5000 %Altitude (ft)
    M = 0.201 %Mach Number
    U_1 = 220.1 %TAS (ft/s)
    q_bar = 49.6 %Dynamic Pressure (lbs/ft^2)
    c_frac = 26.4 %CG Location
    a_1 = 0 %Angle of Attack (Deg)

    %%%%%%%%%%%%%%%%Mass Data For All Flight Conditions%%%%%%%%%%%%%%%%%%%%%%
    %Mass Data in Climb, Cruise, and Approach Conditions
    W = 2650 %Weight (lbs)
    m = W/g %Mass
    I_xx_B = 948 %Moment of Inertia along the X-axis (slug-ft^2)
    I_yy_B = 1346 %Moment of Inertia along the Y-axis (slug-ft^2)
    I_zz_B = 1967 %Moment of Inertia along the Z-axis (slug-ft^2)
    I_xz_B = 0 %Moment of Inertia along the XZ plane (slug-ft^2)

    %%%%%%%%%%%%%%%%%%%%%%%Longitudinal Control and Hinge Momoment $$$$$$$$$

    %Cruise Longitudinal Control and Hinge Moment Derivatives
    Ch_a = -0.0584
    Ch_del_e = -0.585

    %%%%%%%%%%%%%%%%%%%Lateral-Directional Stability%%%%%%%%%%%%%%%%%%%%

    %Cruise Lateral-Directional Stability Derivatives
    Cl_B = -0.0923
    Cl_p = -0.484
    Cl_r = 0.0798
    Cy_B = -0.393
    Cy_p = -0.075
    Cy_r = 0.214
    Cn_B = 0.0587
    Cn_TB = 0
    Cn_p = -0.0278
    Cn_r = -0.0937

    %%%%%%%%%%%%%%%%%%%%Lateral-Directional Control and Hing Moment%%%%%%%%%

    %Cruise Lateral-Directional Control and Hinge Moment Derivatives
    Cl_dela = 0.229
    Cl_delr = 0.0147
    Cy_dela = 0
    Cy_delr = 0.187
    Cn_dela = -0.0216
    Cn_delr = -0.0645
    Ch_aa = 0
    Ch_dela = -0.363
    Ch_Br = 0.0819
    Ch_delr = -0.567

    %Lateral-Directional Dimensional Stability Derivatives
    YB = (q_bar*S*Cy_B)/m
    Yp = (S*b*Cy_p)/(2*m*U_1)
    Yr = (q_bar*S*b*Cy_r)/(2*m*U_1)
    Y_dela = (q_bar*S*Cy_dela)/m
    Y_delr = (q_bar*S*Cy_delr)/m
    LB = (q_bar*S*b*Cl_B)/I_xx_B
    Lp = (q_bar*S*(b^2)*Cl_p)/(2*I_xx_B*U_1)
    Lr = (q_bar*S*(b^2)*Cl_r)/(2*I_xx_B*U_1)
    L_dela = (q_bar*S*b*Cl_dela)/I_xx_B
    L_delr = (q_bar*S*b*Cl_delr)/I_xx_B
    NB = (q_bar*S*b*Cn_B)/I_zz_B
    NT_B = (q_bar*S*b*Cn_TB)/I_zz_B
    Np = (q_bar*S*(b^2)*Cn_p)/(2*I_zz_B*U_1)
    Nr = (q_bar*S*(b^2)*Cn_r)/(2*I_zz_B*U_1)
    N_dela = (q_bar*S*b*Cn_dela)/I_zz_B
    N_delr = (q_bar*S*b*Cn_delr)/I_zz_B

    %Perturbed Lateral-Directional Equations of Motion
    A_1 = I_xz/I_xx_B
    B_1 = I_xz/I_zz_B

    %Denominator Polynomial
    A2 = U_1(1-A_1*B_1)
    B2 = -YB*(1-A_1*B_1)-U_1*(Lp+Nr+A_1*Np+B_1*Lr)
    C2 = U_1*(Lp*Nr-Lr*Np)+YB*(Nr+Lp+A_1*Np+B_1*Lr)-Yp*(LB+NB*A_1+NT_B*A_1)...
        +U_1*(LB*B_1+NB+NT_B)-Yr*(LB*B_1+NB+NT_B)
    D2 = -YB*(Lp*Nr-Lr*Np)+Yp*(LB*Nr-NB*Lr-NT_B*Lr)-g*cos(theta1)...
        *(LB+NB*A_1+NT_B*A_1)+U_1*(LB*Np-NB*Lp-NT_B*Lp)-Yr*(LB*Np-NB*Lp-NT_B*Lp)
    E2 = g*cos(theta1)*(LB*Nr-NB*Lr-NT_B*Lr)

    D_2 = [A2 B2 C2 D2 E2]

    %Numerator Polynomial
    AB = Y_dela*(1-A_1*B_1)+Y_delr*(1-A_1*B_1)
    BB = -Y_dela*(Nr+Lp+A_1*Np+B_1*Lr)-Y_delr*(Nr+Lp+A_1*Np+B_1*Lr)+Yp*(L_dela+N_dela*A_1)+Yp*(L_delr+N_delr*A_1)...
        +Yr*(L_dela*B_1+N_dela)+Yr*(L_delr*B_1+N_delr)-U_1*(L_dela*B_1+N_dela)...
        +U_1*(L_delr*B_1+N_delr)
    CB = Y_dela*(Lp*Nr-Np*Lr)+Y_delr*(Lp*Nr-Np*Lr)+Yp*(N_dela*Lr-L_dela*Nr)+Yp*(N_delr*Lr-L_delr*Nr)...
        +g*cos(theta1)*(L_dela+N_dela*A_1)+g*cos(theta1)*(L_delr+N_delr*A_1)...
        +Yr*(L_dela*Np-N_dela*Lp)+Yr*(L_delr*Np-N_delr*Lp)+U_1*(L_dela*Np-N_dela*Lp)+U_1*(L_delr*Np-N_delr*Lp)
    DB = g*cos(theta1)*(N_dela*Lr-L_dela*Nr)+g*cos(theta1)*(N_delr*Lr-L_delr*Nr)
    N_B = [AB BB CB DB]

    %Transfer Function
    cruise_sys = tf(N_B,D_2)
    damp(cruise_sys)
end
function approach_sys
    g = 32.174
    theta1 = 0
    %Cessna design parameters
    S = 174 %Planform Area (ft^2)
    c = 4.9 %Chord length (ft)
    b = 36.0 %Wing span (ft)

    %%%%%%%%%%%%%%%%Flight Conditions%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %Approach Flight Condition
    h = 0
    M = 0.096
    U_1 = 107.1
    q_bar = 13.6
    c_frac = 26.4
    a_1 = 4

    %%%%%%%%%%%%%%%%Mass Data For All Flight Conditions%%%%%%%%%%%%%%%%%%%%%%
    %Mass Data in Climb, Cruise, and Approach Conditions
    W = 2650 %Weight (lbs)
    m = W/g %Mass
    I_xx_B = 948 %Moment of Inertia along the X-axis (slug-ft^2)
    I_yy_B = 1346 %Moment of Inertia along the Y-axis (slug-ft^2)
    I_zz_B = 1967 %Moment of Inertia along the Z-axis (slug-ft^2)
    I_xz_B = 0 %Moment of Inertia along the XZ plane (slug-ft^2)

    %%%%%%%%%%%%%%%%%%%%%%%Longitudinal Control and Hinge Momoment $$$$$$$$$

    %Approach Longitudinal Control and Hinge Moment Derivatives
    Ch_a = -0.0549
    Ch_del_e = -0.594

    %%%%%%%%%%%%%%%%%%%Lateral-Directional Stability%%%%%%%%%%%%%%%%%%%%

    %Approach Lateral-Directional Stability Derivatives
    Cl_B = -0.0969
    Cl_p = -0.494
    Cl_r = 0.2039
    Cy_B = -0.303
    Cy_p = -0.213
    Cy_r = 0.201
    Cn_B = 0.0701
    Cn_TB = 0
    Cn_p = -0.0960
    Cn_r = -0.1151

    %%%%%%%%%%%%%%%%%%%%Lateral-Directional Control and Hing Moment%%%%%%%%%

    %Approach Lateral-Directional Control and Hinge moment Derivatives
    Cl_dela = 0.229
    Cl_delr = 0.0147
    Cy_dela = 0
    Cy_delr = 0.187
    Cn_dela = -0.0504
    Cn_delr = -0.0805
    Ch_aa = 0
    Ch_dela = -0.369
    Ch_Br = 0.0819
    Ch_delr = -0.579

    %Lateral-Directional Dimensional Stability Derivatives
    YB = (q_bar*S*Cy_B)/m
    Yp = (S*b*Cy_p)/(2*m*U_1)
    Yr = (q_bar*S*b*Cy_r)/(2*m*U_1)
    Y_dela = (q_bar*S*Cy_dela)/m
    Y_delr = (q_bar*S*Cy_delr)/m
    LB = (q_bar*S*b*Cl_B)/I_xx_B
    Lp = (q_bar*S*(b^2)*Cl_p)/(2*I_xx_B*U_1)
    Lr = (q_bar*S*(b^2)*Cl_r)/(2*I_xx_B*U_1)
    L_dela = (q_bar*S*b*Cl_dela)/I_xx_B
    L_delr = (q_bar*S*b*Cl_delr)/I_xx_B
    NB = (q_bar*S*b*Cn_B)/I_zz_B
    NT_B = (q_bar*S*b*Cn_TB)/I_zz_B
    Np = (q_bar*S*(b^2)*Cn_p)/(2*I_zz_B*U_1)
    Nr = (q_bar*S*(b^2)*Cn_r)/(2*I_zz_B*U_1)
    N_dela = (q_bar*S*b*Cn_dela)/I_zz_B
    N_delr = (q_bar*S*b*Cn_delr)/I_zz_B

    %Perturbed Lateral-Directional Equations of Motion
    A_1 = I_xz/I_xx_B
    B_1 = I_xz/I_zz_B

    %Denominator Polynomial
    A2 = U_1(1-A_1*B_1)
    B2 = -YB*(1-A_1*B_1)-U_1*(Lp+Nr+A_1*Np+B_1*Lr)
    C2 = U_1*(Lp*Nr-Lr*Np)+YB*(Nr+Lp+A_1*Np+B_1*Lr)-Yp*(LB+NB*A_1+NT_B*A_1)...
        +U_1*(LB*B_1+NB+NT_B)-Yr*(LB*B_1+NB+NT_B)
    D2 = -YB*(Lp*Nr-Lr*Np)+Yp*(LB*Nr-NB*Lr-NT_B*Lr)-g*cos(theta1)...
        *(LB+NB*A_1+NT_B*A_1)+U_1*(LB*Np-NB*Lp-NT_B*Lp)-Yr*(LB*Np-NB*Lp-NT_B*Lp)
    E2 = g*cos(theta1)*(LB*Nr-NB*Lr-NT_B*Lr)

    D_2 = [A2 B2 C2 D2 E2]

    %Numerator Polynomial
    AB = Y_dela*(1-A_1*B_1)+Y_delr*(1-A_1*B_1)
    BB = -Y_dela*(Nr+Lp+A_1*Np+B_1*Lr)-Y_delr*(Nr+Lp+A_1*Np+B_1*Lr)+Yp*(L_dela+N_dela*A_1)+Yp*(L_delr+N_delr*A_1)...
        +Yr*(L_dela*B_1+N_dela)+Yr*(L_delr*B_1+N_delr)-U_1*(L_dela*B_1+N_dela)...
        +U_1*(L_delr*B_1+N_delr)
    CB = Y_dela*(Lp*Nr-Np*Lr)+Y_delr*(Lp*Nr-Np*Lr)+Yp*(N_dela*Lr-L_dela*Nr)+Yp*(N_delr*Lr-L_delr*Nr)...
        +g*cos(theta1)*(L_dela+N_dela*A_1)+g*cos(theta1)*(L_delr+N_delr*A_1)...
        +Yr*(L_dela*Np-N_dela*Lp)+Yr*(L_delr*Np-N_delr*Lp)+U_1*(L_dela*Np-N_dela*Lp)+U_1*(L_delr*Np-N_delr*Lp)
    DB = g*cos(theta1)*(N_dela*Lr-L_dela*Nr)+g*cos(theta1)*(N_delr*Lr-L_delr*Nr)
    N_B = [AB BB CB DB]

    %Transfer Function
    appoach_sys = tf(N_B,D_2)
    damp(appoach_sys)
end

[Sam Chak]

Hi @Jeffrey Lewis

I actually added the step function in the code to generate step response for each flight condition. Step response always start

https://www.mathworks.com/help/ident/ref/lti.step.html

You can also try the lsim(sys, u, t, x0) function for a more realistic simulated response data because you can design the the input signal u, and specify the initial condition x0. However, you must convert sys from transfer function into a state-space model (preferably in controllable canonical form)

https://www.mathworks.com/help/ident/ref/lti.lsim.html