Error: System does not seem to contain a first-order differential equation for the field when using odeToVectorField

Hello,

I have been trying to work around a way to get this system of ODEs through odeToVectorField, but it keeps returning the following error message:

Error using mupadengine/feval_internal

System does not seem to contain a first-order differential equation for the field 'D(r)(t)'.

If it is any help, the system of ODEs I am trying to solve is the geodesic equations of a particle within the Schwarszchild solution (General Relativity).

The code used is the following, although there is a long and tedious part after defining the symbolic variables I omit, whose purpose is to build the terms used in the ODEs. To make things easier, I proceed to paste the values of x and gds which are used to build the ODEs; and can be checked to be correct with the theoretically obtained geodesic equations, so there should be no problem with them. Also, numeric substitution is not necessary, at least not for using odeToVectorField.

clear all
close
clc
syms ct(s) c t(s) r(s) theta(s) phi(s) M G rs  
%[Here should be the construction of the gds term which is long and complicated, so I just paste the result. x is defined as it is written although some irrelevant substitutions must be done in order to obtain gds so I just paste the final form it will take.]

gds = [-(c*rs*diff(r(s), s)*diff(t(s), s))/(r(s)*(rs - r(s))), (rs - r(s))*diff(theta(s), s)^2 + sin(theta(s))^2*(rs - r(s))*diff(phi(s), s)^2 + (rs*diff(r(s), s)^2)/(2*r(s)*(rs - r(s))) - (c^2*rs*(rs - r(s))*diff(t(s), s)^2)/(2*r(s)^3), (2*diff(r(s), s)*diff(theta(s), s))/r(s) - (sin(2*theta(s))*diff(phi(s), s)^2)/2, (2*diff(phi(s), s)*diff(r(s), s))/r(s) + (2*cos(theta(s))*diff(phi(s), s)*diff(theta(s), s))/sin(theta(s))];
x = [c*t(s),r(s),theta(s),phi(s)];

eqn = [diff(x(1),s,2) == -gds(1),diff(x(2),s,2) == -gds(2), ...
    diff(x(3),s,2) == -gds(3),diff(x(4),s,2) == -gds(4)];

cn = 299792458;
Gn = 6.67408e-11;
Mn = 1.5988e30;
rsn = 2*Gn*Mn/cn^2;
RM = 6.982e10;
mv = 38.86e3;
eqn=subs(eqn,{c,rs},{cn,rsn});
V = odeToVectorField(eqn);

Thank you in advance for taking your time. I must say this is my first time asking, so feel free to make any criticism if necessary.

[Mario Misas Arcos]

SOLVED The problem was that the input I was using for odeToVectorField was a symbolic expression when what it requiered was a symbolic function. Hope this helps anyone.