function [y]=fun(x,u,parameter)
% force balance in the vdot and hdot directions
%   for the MPX5, sea level, weight= 19.2 lbf
% y(1)=mass*vdot
% y(2)=mass*hdotdot
% x(1)=v(ft/sec)
% x(2)=hdot(ft/sec)
% u(1)=alpha(rad)
% u(2)=Bhp(hp)

[n,m]=size(x);
y=zeros(n,1); 

% states
% x(1)=v(ft/sec)
% x(2)=hdot(ft/sec)
v=x(1);
hdot=x(2);
singamma=hdot/v;
gamma=asin(singamma);

% controls
% u(1)=alpha(rad)
% u(2)=Bhp(hp)
alpha=u(1);
Bhp=u(2);


% parameters
cl0=parameter(1);
cla=parameter(2);
cd0=parameter(3);
k=parameter(4);
weight=parameter(5); % lbf
s=parameter(6);     % ft^2
rho=parameter(7); % for sea level
phit=parameter(8);
EtaP=parameter(9);

% nondimensional aerodynamics
cl=cl0 + cla*alpha;
cd=cd0 + k*cl*cl;

%dimensional aerodynamics
qs=.5*rho*v*v*s;
lift=qs*cl;
drag=qs*cd;
thrust=550*Bhp*EtaP/v;

cosgamma=cos(gamma);
y(1)=thrust*cos(phit+alpha) -drag - weight*singamma;
y(2)=lift*cosgamma + thrust*sin(phit+alpha+gamma) -weight -drag*singamma;
% y(1)=mass*vdot
% y(2)=mass*hdotdot
% x(1)=v(ft/sec)
% x(2)=hdot(ft/sec)
% u(1)=alpha(rad)
% u(2)=Bhp(hp)

