Dont bother with the hassle to go through symbolic-variables and ode2vectorfield when transforming your 2 clean 2nd-order ODEs. Convert them by hand!
function dy12dt = myODE(t,y,mu)
y1 = y(1);
y2 = y(2);
dy1dt = y(3);
dy2dt = y(4);
D1 = ((y1+mu)^2 + y2^2)^(3/2);
d2y1dt2 = y1 + 2*dy2dt + (1-mu)*(y1+mu)/D1;
dy12dt = [dy1dt;dy2dt;d2y1dt2;d2y2dt2]
Then you have a system of four coupled first-order ODEs to solve with ode45 with time-span and initial conditions of your choice...