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
