Difficulty simulating trajectory of moon and sun for Apollo 11 mission

Question

As an engineering student, I am doing a project about Kalman filtering and the mission Apollo 11. Now I tried to simulate the trajectory of the moon and sun for 3 days (the duration in which the spacecraft flight from the earth to the moon). But to my surprise, when I solved the motion equations of the spacecraft, I did not get good results.

How does this trajectory look like?

Here are blocks of the code:

function dydt = ode(t, y)
% Extract the state variables
X = y(1);
Y = y(2);
Z = y(3);
Xdot = y(4);
Ydot = y(5);
Zdot = y(6);
%% date du départ 16 juillet 1969
mois = 7;
jours = [16 17 18];
ann = 1969;
%d = 367Y - (7(Y + ((M+9)/12)))/4 + (275*M)/9 + D - 730530
d1 = 367ann - (7(ann + ((mois+9)/12)))/mois + (275*mois)/9 + jours(1) - 
730530;
d2 = 367ann - (7(ann + ((mois+9)/12)))/mois + (275*mois)/9 + jours(2) - 730530;
d3 = 367ann - (7(ann + ((mois+9)/12)))/mois + (275*mois)/9 + jours(3) - 730530;
%% sun position
if t <= 2460
    d = d1;
elseif t>2460 & t <= 4860
    d = d1 + (t - 2460) / (4860 - 2460) * (d2 - d1);
else
    d = d2 + (t - 4860) / (7260 - 4860)  (d3 - d2);
end
w = 282.9404 + 4.70935e-5* d ;%   (longitude of perihelion)
a = 1.000000    ;   %   (mean distance, a.u.)
e = 0.016709 - 1.15e-9* d ; %%  (eccentricity)
M = 356.0470 + 0.9856002585* d ; %  (mean anomaly)
oblecl = 23.4393 - 3.563e-7 * d;
L = w + M ;
E = M + (180/pi) * e * sin(M) * (1 + e * cos(M));
x = cos(E) - e;
y = sin(E) * sqrt(1 - e*e);
r = sqrt(xx + yy);
v = atan2( y, x );
lon = v + w;
x = r * cos(lon);
y = r * sin(lon);
z = 0.0;
xequat = 0.881048;
yequat = 0.482098 * cos(oblecl) - 0.0 * sin(oblecl);
zequat = 0.482098 * sin(oblecl) + 0.0 * cos(oblecl);
Xs = xequat*1.496e+11;
Ys = yequat*1.496e+11;
Zs = zequat*1.496e+11;
%% lune position     !!!!
N = 125.1228 - 0.0529538083 * d;   % (Long asc. node)
inc =   5.1454;                          % (Inclination)
w = 318.0634 + 0.1643573223 * d;   % (Arg. of perigee)
a =  60.2666  ;                             % (Mean distance)
e = 0.054900  ;                             % (Eccentricity)
M = 115.3654 + 13.0649929509 * d;    %(Mean anomaly)
% normalisation
M = M + 129*360;
w = w + 360;
E0 = M + (180/pi) * e * sin(M) * (1 + e * cos(M));
E1 = E0 - (E0 - (180/pi) * e * sin(E0) - M) / (1 - e * cos(E0));
E = E1;
x = a * (cos(E) - e);
y = a * sqrt(1 - ee)  sin(E);
r = sqrt( xx + yy );
v = atan2( y, x );
xeclip = r * ( cos(N) * cos(v+w) - sin(N) * sin(v+w) * cos(inc) );
yeclip = r * ( sin(N) * cos(v+w) + cos(N) * sin(v+w) * cos(inc) );
zeclip = r * sin(v+w) * sin(inc);
xeq = xeclip;
yeq = yeclip * cos(oblecl) - zeclip * sin(oblecl);
zeq = yeclip * sin(oblecl) + zeclip * cos(oblecl);
Xm = xeq*1.496e+11 ;
Ym = yeq*1.496e+11;
Zm = zeq*1.496e+11;
%% constantes
mu_e = 3.986135e14; %sets the gravitational parameter for the Earth.
mu_m = 4.89820e12;%sets the gravitational parameter for the Moon.
mu_s = 1.3253e20;%sets the gravitational parameter for the Sun.
a = 6.37826e6; %the equatorial radius of the Earth.
J = 1.6246e-3;%sets the Earth's second dynamic form factor.
% Calculate the distances and other parameters
r = sqrt(X^2 + Y^2 + Z^2);
rs = sqrt(Xs^2 + Ys^2 + Zs^2);
rm = sqrt(Xm^2 + Ym^2 + Zm^2);
delta_m = sqrt((X - Xm)^2 + (Y - Ym)^2 + (Z - Zm)^2);
delta_s = sqrt((X - Xs)^2 + (Y - Ys)^2 + (Z - Zs)^2);
Y_s = Ys;
Z_s = Zs;
% Calculate the derivatives
X2dot = -((mu_eY)/r^3)  (1 + (J(a/r)^2)(1 - 5*Z^2/r^2)) ...
- (mu_m*(X - Xm)/delta_m^3) ...

- (mu_m*Xm/rm^3) ...

- (mu_s*(X - Xs)/delta_s^3) ...

- (mu_s*Xs/rs^3);


Y2dot = -((mu_eY)/r^3)  (1 + (J(a/r)^2)(1 - 5*Z^2/r^2)) ...
- (mu_m*(Y - Ym)/delta_m^3) ...

- (mu_m*Ym/rm^3) ...

- (mu_s*(Y - Y_s)/delta_s^3) ...

- (mu_s*Ys/rs^3);


Z2dot = -((mu_eZ)/r^3)  (1 + (J(a/r)^2)(3 - 5*Z^2/r^2)) ...
- (mu_m*(Z - Zm)/delta_m^3) ...

- (mu_m*Zm/rm^3) ...

- (mu_s*(Z - Z_s)/delta_s^3) ...

- (mu_s*Zs/rs^3);



% Return the derivative vector
%X2dot = interp1(Xdot,X2dot,t);
dydt = [Xdot; Ydot; Zdot; X2dot; Y2dot; Z2dot];
end
Then we have the ODE45 to resolve the equation
close all
clear all; clc;
%% Set the initial conditions
X0 =100 + 6.37826e6;

Y0 = 100 + 6.37826e6;

Z0 = 100 + 6.37826e6; 

 Xdot0 =  -2.0225e+14;    %Xm/72*3600, repère géocentrique %8.9919e+10  ;

 Ydot0 =  -1.8023e+14;%-1.0908e+11;

 Zdot0 =  3.4945e+14; %-4.8960e+10;

y0 = [X0; Y0; Z0; Xdot0; Ydot0; Zdot0];

tspan = 0:1:72*60*60;


[t, y] = ode45(@ode, tspan, y0);
plot3(y(:,1),y(:,2),y(:,3))
title('Solution of the Motion Equation with ODE45');