Error using sym function with dsolve

So I am trying to implement a 3DOF model to plot reentry trajectories but I keep getting the following error.

Error using sym/subsindex (line 769)
Invalid indexing or function definition. When defining a function, ensure that the arguments are
symbolic variables and the body of the function is a SYM expression. When indexing, the input
must be numeric, logical, or ':'.

My code is as follows:

syms r(t) theta(t) phi(t) V(t) gamma(t) chi(t)
%r,theta,phi=radius,longitude,latitude
%gamma=flight path angle
%chi = flight path azimute
%V=velocity

D=1 %Drag force (N)
m=1 %mass (kg)
mu=1 %graviational constant
R=1%????
J=1%second gravitational harmonics
RE=1%????
gr=-mu/R^2*(1-3/2*J*(RE/R)^2*(3*(cos(phi))^2-1)); %radial component of acceleration due to gravite
gphi=3*mu/R^2*(RE/R)^2*sin(phi)*cos(phi)*J; %latitudinal compenent of g
sigma=1 %
L= 1%lift

%%dsolve

ode1= diff(r)== -sin(gamma)*V;
ode2= diff(theta)==cos(gamma)*cos(chi)/(r*cos(phi))*V;
ode3= diff(phi) == cos(gamma)*sin(chi)/r*V;
ode4= diff(V) == -D/m+gr*sin(gamma)-gphi*cos(gamma)*cos(chi);
ode5=diff(gamma) == 1/V*(gr*cos(gamma)+gphi*sin(gamma)*cos(chi)-V^2/r*cos(gamma)-L/m*cos(sigma));
ode6=diff(chi)==-cos(gamma)*cos(chi)*tan(phi)/r*V+gphi*cos(gamma)*sin(chi)/V+L/m*cos(gamma)*sin(sigma)/V;

odes=[ode1;ode2;ode3;ode4;ode5;ode6];

[rsol(t),thetasol(t),phisol(t),Vsol(t),gammasol(t),chisol(t)]=dsolve(odes)

Answer 1

[Walter Roberson]

[rsol,thetasol,phisol,Vsol,gammasol,chisol]=dsolve(odes)

However, the system is too complicated for MATLAB. Other packages struggle with it as well.