how to include variables for plotting

Hi, I have this function and script files. I would like to know how I could plot the control moments?? de, da, dr not sure how to do it. Thanks, any help will be gretly appreciated.

function xdot = STOL_EOM(t,x)
% STOL_EOM contains the nonlinear equations of motion for a rigid airplane.
% (NOTE: The aerodynamic model is linear.)
global e1 e2 e3 rho m g S b c AR WE Power ...
    GeneralizedInertia GeneralizedInertia_Inv ...
    CD0 e CL_Trim CLalpha CLq CLde Cmde ...
    CYbeta CYp CYr Cmalpha Cmq Clbeta Clp Clr Cnbeta Cnp Cnr ...
    CYda CYdr Clda Cldr Cnda Cndr ...
    de0 da0 dr0 
X = x(1:3);
Theta = x(4:6);
    phi = Theta(1);
    theta = Theta(2);
    psi = Theta(3);
V = x(7:9);
    u = V(1);
    v = V(2);
    w = V(3);    
omega = x(10:12);
    p = omega(1);
    q = omega(2);
    r = omega(3);    
alpha = atan(w/u);
beta = asin(v/norm(V));
P_dynamic = (1/2)*rho*norm(V)^2;
RIB = expm(psi*hat(e3))*expm(theta*hat(e2))*expm(phi*hat(e1));
LIB = [1, sin(phi)*tan(theta), cos(phi)*tan(theta);
    0, cos(phi), -sin(phi);
    0, sin(phi)/cos(theta), cos(phi)/cos(theta)];
VE = V + RIB'*WE;
% Kinematic equations
XDot = RIB*VE;
ThetaDot = LIB*omega;
% Weight
W = RIB'*(m*g*e3);
% Control moments
de = de0;
da = da0;
dr = dr0;
% Components of aerodynamic force (modulo "unsteady" terms)
CL = CL_Trim + CLalpha*alpha + CLq*((q*c)/(2*norm(V))) + CLde*de;
CD = CD0 + (CL^2)/(e*pi*AR);
temp = expm(-alpha*hat(e2))*expm(beta*hat(e3))*[-CD; 0; -CL];
CX = temp(1);
CZ = temp(3);
CY = CYbeta*beta + CYp*((b*p)/(2*norm(V))) + CYr*((b*r)/(2*norm(V))) + ...
    CYda*da + CYdr*dr;
%X = P_dynamic*S*CX + T0;  %Constant Thrust?
Thrust = Power/norm(V);
X = P_dynamic*S*CX + Thrust;  %Constant Thrust?
Y = P_dynamic*S*CY;
Z = P_dynamic*S*CZ;
Force_Aero = [X; Y; Z];
% Components of aerodynamic moment (modulo "unsteady" terms)
Cl = Clbeta*beta + Clp*((b*p)/(2*norm(V))) + Clr*((b*r)/(2*norm(V))) + ...
    Clda*da + Cldr*dr;
Cm = Cmalpha*alpha + Cmq*((c*q)/(2*norm(V))) + Cmde*de;
Cn = Cnbeta*beta + Cnp*((b*p)/(2*norm(V))) + Cnr*((b*r)/(2*norm(V))) + ...
    Cnda*da + Cndr*dr;
L = P_dynamic*S*b*Cl;
M = P_dynamic*S*c*Cm;
N = P_dynamic*S*b*Cn;
Moment_Aero = [L; M; N];
% Sum of forces and moments
Force = W + Force_Aero;
Moment = Moment_Aero;
% NOTE: GeneralizedInertia includes "unsteady" aerodynamic terms
% (i.e., added mass/inertia).
temp = GeneralizedInertia*[VE; omega];
LinearMomentum = temp(1:3);
AngularMomentum = temp(4:6);
clear temp
RHS = [cross(LinearMomentum,omega) + Force; ...
    cross(AngularMomentum,omega) + Moment];
temp = GeneralizedInertia_Inv*RHS;
VEDot =  temp(1:3);
VDot = VEDot + cross(omega,RIB'*WE);
omegaDot = temp(4:6);
clear temp
xdot = [XDot; ThetaDot; VDot; omegaDot];
====================================script=================
clear
close all
% GenericFixedWingScript.m solves the nonlinear equations of motion for a 
% fixed-wing aircraft flying in ambient wind with turbulence. 
global e1 e2 e3 rho m g S b c AR Inertia 
% Basis vectors
e1 = [1;0;0];
e2 = [0;1;0];
e3 = [0;0;1];
% Atmospheric and gravity parameters (Constant altitude: Sea level)
%a = 340.3;          % Speed of sound (m/s)
rho = 1.225;        % Density (kg/m^3)
g = 9.80665;        % Gravitational acceleration (m/s^2)
%u0 = Ma*a;
u0 = 13.7;
P_dynamic = (1/2)*rho*u0^2;
% Aircraft parameters (Bix3)
m = 1.202;          % Mass (slugs)
W = m*g;            % Weight (Newtons)
Ix = 0.2163;         % Roll inertia (kg-m^2)
Iy = 0.1823;         % Pitch inertia (kg-m^2)
Iz = 0.3396;         % Yaw inertia (kg-m^2)
Ixz = 0.0364;            % Roll/Yaw product of inertia (kg-m^2)
%Inertia = [Ix, 0, -Ixz; 0, Iy, 0; -Ixz, 0, Iz];
S = 0.0285;         % Wing area (m^2)
b = 1.54;           % Wing span (m)
c = 0.188;          % Wing chord (m)
AR = (b^2)/S;       % Aspect ratio
CL_Trim = W/(P_dynamic*S);
CD_Trim = 0.036;
e = 0.6;    %Contrived
CD0 = CD_Trim - (CL_Trim^2)/(e*pi*AR);
% Equilibrium power (constant)
Thrust_Trim = P_dynamic*S*CD_Trim;
Power = Thrust_Trim*u0;
% Longitudinal nondimensional stability and control derivatives
%Cx0 = 0.197;
%Cxu = -0.156;
%Cxw = 0.297;
%Cxw2 = 0.960;
%Cz0 = -0.179;
%Czw = -5.32;
%Czw2 = 7.02;
%Czq = -8.20;
%Czde = -0.308;
%Cm0 = 0.0134;
%Cmw = -0.240;
%Cmq = -4.49;
%Cmde = -0.364;
% Lateral-directional nondimensional stability and control derivatives
X0 = ones(3,1);
Theta0 = ones(3,1);
V0 = u0*e1;
omega0 = ones(3,1);
y0 = [X0; Theta0; V0; omega0];
t_final = 10;
[t,y] = ode45('STOL_EOM',[0:0.1:t_final]',y0);
figure(1)
subplot(2,1,1)
plot(t,y(:,1:2))
ylabel('Position')
subplot(2,1,2)
plot(t,y(:,3:4))
ylabel('Attitude')
figure(2)
subplot(2,1,1)
plot(t,y(:,6:8))
ylabel('Velocity')
subplot(2,1,2)
plot(t,y(:,9:11))
ylabel('Angular Velocity')
Cesar Cardenas