How to solve system of 2nd order differential equations using ode45

I have three 2nd order differential equations with my initial conditions and I'm trying to use the ode45 function in matlab to solve this.

I wish get (x,y) position in x-y plane

but i can`t simultaneously get x-position and y-position respect t.

%planar circular restricted 3 body problem

% %Initial Conditions

% x(0) = 35.0; %AU

% y(0) = 0; %AU

%Constants

G = 1;

m1 = 0.52; %solar mass

m2 = 0.33; %solar mass

r1= -21.3; %AU

r2= 33.7; %AU

ome=sqrt((m1+m2)/(-r1+r2)^3);

% x=x(1) x'=x(2) y=x(3) y'=x(4)

% x''=x*ome^2+2*ome*y'-m1*(x-r1)/(((x-r1)^2+y^2)^1.5)-m2*(x-r2)/(((x-r2)^2+y^2)^1.5)

% y''=y*ome^2-2*ome*x'-m1*y/(((x-r1)^2+y^2)^1.5)-m2*y/(((x-r2)^2+y^2)^1.5)

% x=x1 , x1'=x2 . y=x3 , x3'=x4

[t,x]=ode45(@crtbp,[0 35],[35;0;0;0.03]);

function dxdt=crtbp(t,x)

m1=0.52;

m2=0.33;

r1=-21.3;

r2=33.7;

ome=sqrt((m1+m2)/(-r1+r2)^3);

dxdt=[x(2); x(1)*ome^2+2*ome*x(4)-m1*(x(1)-r1)/(((x(1)-r1)^2+x(3)^2)^1.5)-m2*(x(1)-r2)/(((x(1)-r2)^2+x(3)^2)^1.5); x(4); x(3)*ome^2-2*ome*x(2)-m1*x(3)/(((x(1)-r1)^2+x(3)^2)^1.5)-m2*x(3)/(((x(1)-r2)^2+x(3)^2)^1.5)];

end

[darova]

Looks correct. What is the problem? You want you plot x vs y?

plot(x(:,1),x(:,3))

Answer 1

[Stephan]

plot3(t,x(:,1),x(:,3))