Solving an ode for newton with vectors

Hello;

I am used to use python for basic programming stuff and was now told to solve an ode in matlab, so I am a bit in need of help.

What I understood so far, I made a code solving a simple ODE, r''=-w²r , just to understand the concept:

function odesolver
w=2; %random value
r0=1;%my values for r(t=0)=r0 and r'(t=0)=v0
v0=0;
[T,X]=ode45(@diff,[0 10],[r0 v0]);
plot(T,X)
    function drdt=diff(t,r) %function for the right side
        drdt=[r(2);-w^2*r(1)];
    end
end

So far so good.

Now I don't know if my problem is one or multiple questions, I am just going to ask it:

I have to solve the 2-body newton ODE with vectors. Meaning:

My ODE's right side looks like this:

GM is just a constant as w was before, but the r is a vector, and so are my starting values.

Throwing that into my code leaves me with:

function odesolver
GM=6; %just random numbers so far to get it to work at all
r0=[2;3;3]; %here are now vectors
v0=[3;5;6];
[T,X]=ode45(@right,[0 10],[r0 v0]);
plot(T,X) % I know "plot(...) is only for 2D, so I have to use meshm(..) for 3D?
    function drdt=right(t,r)
        drdt=[r(2);-(GM/(norm(r(1))^3))*r(1)];%here I tried to implement the new right side
    end
end

Obviously that leaves me with some errors, mainly, and by looking through other asked questions it has to do something with the dimensions, that I use two vectors with a total of 6 entries but somehow "ask" the program to only give out two values?

Error using odearguments
UNTITLED3/RIGHT returns a vector of length 2, but the length of initial conditions vector is 6. The vector returned by
UNTITLED3/RIGHT and the initial conditions vector must have the same number of elements.

Error in ode45 (line 107)
  odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Error in untitled3 (line 5)
[T,X]=ode45(@right,[0 10],[r0 v0]);

So my question comes down to: How can I solve an ODE with vectors? Like in my case?

Thanks in advance!

Answer 1

[Bjorn Gustavsson]

The thing you have to fix here is the indexing of r in your equation-of-motion. You have to extract the 3 velcity-components for dr/dt, and use the three spatial coordinates to calculate dv/dt:

function odesolver
GM=6; %just random numbers so far to get it to work at all
r0=[2;3;3]; %here are now vectors
v0=[3;5;6];
[T,X]=ode45(@right,[0 10],[r0;v0]); % Edit! Just noticed this, the initial condition should be 1-by-6 or 6-by-1
plot(T,X) % I know "plot(...) is only for 2D, so I have to use meshm(..) for 3D?
    function drdt=right(t,rv) % changed variable-name to the more explicit RV
        r = rv(1:3); % position-vector is first three components
        v = rv(4:6); % velocity-vector
        drdt=[v(:);-(GM/(norm(r(:))^3))*r(:)];%here I tried to implement the new right side
    end
end

I find it clearer to explicitly name the 6-d r-v vector rv, and then separate the velocty and position-components up top in the equation of motion.

HTH