System of coupled ODEs

Hello,

I'm trying to solve a system of coupled ODEs. The odes are:

x'' + 2y' = d/dx(F(x,y))

y'' - 2x' = d/dy(F(x,y))

where F is some function of x and y and I'm using ode45 solver. My code looks like this:

function Ord2ODE

t=0:0.0001:10;
%x,xdot,y,ydot
ainit = [50; 0; 0; 20];

[t,a] = ode45(@rhs,t,ainit);

figure;
plot(a(:,1), a(:,2));

end

where function 'rhs' is:

function dadt = rhs(t,a)

mu = 0;

x = a(1);
xdot = a(2);
y = a(3);
ydot = a(4);

Fy = %formula depending on a(1), a(2)
Fx = %formula depending on a(1), a(2)

dadt1 = a(2);
dadt3 = a(4);
dadt2 = 2*a(4) - Fx; 
dadt4 = -2*a(2) - Fy;

dadt = [dadt1;dadt2;dadt3;dadt4];

end

There's plenty of examples of solving uncoupled ODEs on the Internet, but will this work for a set of coupled ones? This seems a little bit funny, since my simulation is for certain gravitational interaction, and x,y are coordinates of an orbiting body. If I set my initial conditions so that initial velocities are all zero, I get a trajectory that is a spiral with increasing coil size - which surely isn't correct, as the object should be bounded. I tried ridicolously small steps and staying away from the mass in the centre to avoid problems with precision, but I still get the same outcome. These sanity checks suggest, taht there is somethign wrong about my code. Any ideas? I would be grateful for any sort of help.

Answer 1

[Minh Tran]

just solve them like any other system of first-order ODE,

but first, you need to convert coupled ones to first-order onesones,

looks like you done it right.