3-dimensional 2nd order ODE, for plotting orbit (integrating gravity equation) on MATLAB

Question

I'd like to ask a question about a problem that has been bothering me for a few days.

I want to integrate a very simple gravitational equation to find the X, Y, Z (position) of an object. It's similar to solving a second-order ODE with three axes X, Y, and Z.

The desired result is to have an orbit that repeats like a spring (as shown in the picture). However, the result from my Matlab code is more like a projectile trajectory.

I think all the lines are appropriate. Therefore, I'm confused as to why the accurate solution is not emerging. Is there any issue that I might be overlooking?

I would appreciate it if you could point out the issues in my code. The code I have written can be seen below.

Gravity Equation

correct graph

my wrong graph

clc; clear;

%% parameters
G = 6.67430e-11;  % J*m/kg^2 Gravitational constant
M = 5.9722e24; % kg, Earth Mass

%% Initial conditions
% x0 y0 z0 xdot0 ydot0 zdot0
initial_state = [6.0827e6; 3.4907e6; 4.6517e6; -651.3041; 1082.1; 7426.4];  % Example initial state

%% time
n = 0.001078; % rad/s
p = 2*pi/n; % period
tspan = [0 p];  % Time span for the simulation

%% Solve
[t, result] = ode45(@(t, y) unconstranined_ECI_ode(t, y, G, M), tspan, initial_state);

% ECI X Y Z
X = result(:, 1);
Y = result(:, 2);
Z = result(:, 3);
ECImatrix = [X Y Z];

%% Plot
figure;

subplot(1, 2, 1);
plot(X, Z);
title('Unconstrained Motion in the xz-plane');
xlabel('X');
ylabel('Z');
grid on;

subplot(1, 2, 2)
plot(Y, Z);
title('Unconstrained Motion in the yz-plane');
xlabel('Y');
ylabel('Z');
grid on;

%% ODE solve function
function state = unconstranined_ECI_ode(t, y, G, M)

    vel = y(4:6); % X dot
    r = y(1:3); % X
    rnorm = norm(r);
    rhat = r/rnorm;
    accel = -(G*M/rnorm^2)*rhat;
    state = [vel;accel];
    
end