问用MATLAB求解二阶微分方程组

EN

x¨(t)=-1/2m ρAC_d cos⁡(arctan⁡((y˙(t))/(x˙(t) )))(〖x˙(t)〗^2+ 〖y˙(t)〗^2)
y¨(t)=-1/2m(2mg+ρAC_d sin⁡(arctan⁡((y˙(t))/(x˙(t) )))(〖x˙(t)〗^2+ 〖y˙(t)〗^2)

# With the initial conditions:

x˙(0)=cosθ ∙ V_0

y˙(0)=sinθ ∙ V_0

x(0)=0

y(0)=0

syms x(t) y(t) a b c d u theta
% Equations
% d2x = -a*cos(arctan(dy/dx))*(((dx)^2)+(dy)^2));
% d2y = -b*(c + d*sin(arctan(dy/dx))*(((dx)^2)+(dy)^2));

%Constants
dx=diff(x,t);
dy=diff(y,t);
d2x=diff(x,t,2);
d2y=diff(y,t,2);
a=1;
b=2;
c=3;
d=4;

%Initial Conditions

cond1 = x(0) == 0;
cond2 = y(0) == 0;
cond3 = dx(0) == u*cos(theta);
cond4 = dy(0) == u*sin(theta);

conds = [cond1 cond2 cond3 cond4];

eq1 = -a*cos(atan(dy/dx))*(((dx)^2)+(dy)^2);
eq2 = -b*(c + d*sin(atan(dy/dx))*(((dx)^2)+(dy)^2));

vars = [x(t); y(t)];
V = odeToVectorField([eq1,eq2]);
M = matlabFunction(V,'vars', {'t','Y'});
interval = [0 5];  %time interval    
ySol = ode23(M,interval,conds);

Error using mupadengine/feval (line 187)
System contains a nonlinear equation in 'diff(y(t), t)'. The system must be quasi-linear:
highest derivatives must enter the differential equations linearly.

Error in odeToVectorField>mupadOdeToVectorField (line 189)
T = feval(symengine,'symobj::odeToVectorField',sys,x,stringInput);

Error in odeToVectorField (line 138)
sol = mupadOdeToVectorField(varargin);

Error in velocity_takeoff (line 29)
V = odeToVectorField([eq1,eq2]);