how to obtain plots from a given function

Hi, not sure how to plot xd1, xd2,,,,,xd3,,,any help will be greatly appreciated. Thanks

function xdot = EoMv2(t,x)

	global mUser Ixx Iyy Izz Ixz S b cBar SMI CONHIS u tuHis deluHis ...
        ACtype MODEL TrimHist RUNNING
    
    D2R =   pi/180;
    R2D =   180/pi;
    
%   Terminate if Altitude becomes Negative    
    [value,isterminal,direction] = event(t,x);
    
%	Earth-to-Body-Axis Transformation Matrix
	HEB		=	DCM(x(10),x(11),x(12));
    
%	Atmospheric State
    x(6)    =   min(x(6),0); % Limit x(6) <= 0 m
	[airDens,airPres,temp,soundSpeed]	=	Atmos(-x(6));
    
%	Body-Axis Wind Field

	windb	=	WindField(-x(6),x(10),x(11),x(12));
    
%	Body-Axis Gravity Components
	gb		=	HEB * [0;0;9.80665];

%	Air-Relative Velocity Vector
    x(1)    =   max(x(1),0);    % Limit axial velocity to >= 0 m/s
	Va		=	[x(1);x(2);x(3)] + windb;
	V		=	sqrt(Va'*Va);
	alphar	=	atan(Va(3)/abs(Va(1)));
    
	alpha 	=	R2D * alphar;
    
	betar	= 	asin(Va(2)/V);
	beta	= 	R2D*betar;
	Mach	= 	V/soundSpeed;
	qbar	=	0.5*airDens*V^2;

%	Incremental Flight Control Effects

	if CONHIS >=1 && RUNNING == 1
		[uInc]	=	interp1(tuHis,deluHis,t);
		uInc	=	(uInc)';
		uTotal	=	u + uInc;
	else
		uTotal	=	u;
    end

%	Force and Moment Coefficients
    if MODEL == 'Alph'
	    [CD,CL,CY,Cl,Cm,Cn,mRef,S,Ixx,Iyy,Izz,Ixz,cBar,b,Thrust] = AlphaModel(x,uTotal,alphar,betar,V);
    end
%     if MODEL == 'Mach'
%         [CD,CL,CY,Cl,Cm,Cn,mRef,S,Ixx,Iyy,Izz,Ixz,cBar,b,Thrust] = MachModel(x,uTotal,Mach,alphar,betar,V);
%     end    

	qbarS	=	qbar*S;

	CX	=	-CD*cos(alphar) + CL*sin(alphar);	% Body-axis X coefficient
	CZ	= 	-CD*sin(alphar) - CL*cos(alphar);	% Body-axis Z coefficient

%	State Accelerations
   
	Xb  =	(CX*qbarS + Thrust)/mRef;
	Yb  =	CY*qbarS/mRef;
	Zb  =	CZ*qbarS/mRef;
	Lb  =	Cl*qbarS*b;
	Mb  =	Cm*qbarS*cBar;
	Nb  =	Cn*qbarS*b;
	nz  =	-Zb/9.80665;      % Normal load factor

%	Dynamic Equations
	xd1 = Xb + gb(1) + x(9)*x(2) - x(8)*x(3);
	xd2 = Yb + gb(2) - x(9)*x(1) + x(7)*x(3);
	xd3 = Zb + gb(3) + x(8)*x(1) - x(7)*x(2);
	
	y	=	HEB' * [x(1);x(2);x(3)];
	xd4	=	y(1);
	xd5	=	y(2);
	xd6	=	y(3);

    xd7	= 	(Izz * Lb + Ixz*Nb - (Ixz*(Iyy - Ixx - Izz)*x(7) + ...
			(Ixz^2 + Izz*(Izz - Iyy))*x(9)) * x(8))/(Ixx * Izz - Ixz^2);
	xd8 = 	(Mb - (Ixx - Izz)*x(7)*x(9) - Ixz*(x(7)^2 - x(9)^2))/Iyy;
	xd9 =	(Ixz*Lb + Ixx*Nb + (Ixz*(Iyy - Ixx - Izz)*x(9) + ...
			(Ixz^2 + Ixx*(Ixx - Iyy))*x(7))*x(8))/(Ixx*Izz - Ixz^2);

	cosPitch	=	cos(x(11));
	if abs(cosPitch)	<=	0.00001
		cosPitch	=	0.00001 * sign(cosPitch);
	end
	tanPitch	=	sin(x(11)) / cosPitch;
		
	xd10	=	x(7) + (sin(x(10))*x(8) + cos(x(10))*x(9))*tanPitch;
	xd11	=	cos(x(10))*x(8) - sin(x(10))*x(9);
	xd12	=	(sin(x(10))*x(8) + cos(x(10))*x(9))/cosPitch;
	
	xdot	=	[xd1;xd2;xd3;xd4;xd5;xd6;xd7;xd8;xd9;xd10;xd11;xd12];
end

[Torsten]

Plotting should be done when the ODE integrator returns a solution, not from the ODE function used to get the derivatives.

[Cesar Cardenas]

Thanks for your reply. Yes, thanks I solved it. Now, I need to apply a square doublet with 1 second of duration. I do not really know how to do it. Also, I need to incorporate this dynamic model:

this is my attempt but not sure,? Any help will be greatly appreciated. Thanks

%Dynamic Actuator Model
 d = (e^-(t*s)/t1*s + 1)*dc;

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')