Satellite x ECI(Km) Satellite y Satellite z 5049.164174 3340.361 3048.609 5227.17629 -2122.7-3755.809 -3260.609868 -4044.645 -4358.112 -6362.512242 755.5592 2229.638 1058.797003 4290.472 5140.661 6728.771277 689.2465-429.4266 1247.121325-4050.981-5289.548 -6308.516462 -2051.021-1421.772 -3435.588934 3351.491 4790.884 5119.921803 3176.862 3103.297 5195.104158-2279.925 -3706.763 -3340.400374-3936.052-4397.115 -6359.488878 948.4093 2164.194 1134.262117 4249.289 5158.093 6750.25117 484.5558-354.7835 1182.490033 -4080.364 -5281.146 -6350.914903 -1858.111-1496.216 -3389.686417 3447.197 4755.526 5174.66388 3019.266 3168.457 5169.875543-2428.657-3646.841 -3399.571419 -3831.134-4443.556

I am trying to plot an orbit in MATLAB. There is something wrong with my code because the final values I get are incorrect. The code is shown below. The correct values are in the image.

mu = 3.986*10^5; % Earth's gravitational parameter [km^3/s^2]

% Transforming orbital elements to cartesian coordinate system for LEO
a_1 = 6782.99;
e_1 = 0.000685539;
inc_1 = 51.64;
v_1 = 5;
argp_1 = 30;
raan_1 = 10;

[x_1, y_1, z_1, vx_1, vy_1, vz_1] = kep2cart(a_1, e_1, inc_1, raan_1, ...
                                            argp_1, v_1);

Y_1 = [x_1, y_1, z_1, vx_1, vy_1, vz_1];

% time_span for two revolutions (depends on the orbit)
t1 = [0 (180*60)];
 

% Setting tolerances
options = odeset('RelTol',1e-12,'AbsTol',1e-12);

% Using ODE45 to numerically integrate for LEO
[t_1, state_1] = ode45(@OrbitProp, t1, Y_1, options);

function dYdt = OrbitProp(t, Y)

mu = 3.986*10^5; % Earth's gravitational parameter [km^3/s^2]

% State Vector

x = Y(1); % [km]

y = Y(2); % [km]

z = Y(3); % [km]

vx = Y(4); % [km/s]

vy = Y(5); % [km/s]

vz = Y(6); % [km/s]

norm = sqrt(x^2+y^2+z^2);

% Differential Equations

vdot1 = ((x)/(norm^3))*(-mu); % [km/s^2]

vdot2 = ((y)/(norm^3))*(-mu); % [km/s^2]

vdot3 = ((z)/(norm^3))*(-mu); % [km/s^2]

% time derivative state vector

dYdt = [vx;vy;vz;vdot1;vdot2;vdot3]; % [km/s and km/s^2]

end

function [x, y, z, vx, vy, vz] = kep2cart(a, e, inc, raan, argp, v)

% Gravitational parameter

mu = 398600.4418; % (km^3/s^2)

% Calculate angular momentum

h = sqrt(a * (1-e^2)*mu); % (km^2/s)

% Calculate position and velocity in the periforcal frame

if e ~= 1

p = a*(1 - e^2);

else

p = h^2/mu;

end

r1 = p*cosd(v) / (1 + e*cosd(v));

r2 = p*sind(v) / (1 + e*cosd(v));

r3 = 0;

v1 = -sqrt(mu/p)*sind(v);

v2 = sqrt(mu/p)*(e + cosd(v));

v3 = 0;

r_w = [r1, r2, r3];

v_w = [v1, v2, v3];

% Define the Rotational Matrices

R1 = [cosd(-raan) sind(-raan) 0;

-sind(-raan) cosd(-raan) 0;

0 0 1];

R2 = [1 0 0;

0 cosd(-inc) sind(-inc);

0 -sind(-inc) cosd(-inc)];

R3 = [cosd(-argp) sind(-argp) 0;

-sind(-argp) cosd(-argp) 0;

0 0 1];

% Calculate position vector

r_rot = (R1 * R2 * R3) * r_w';

% Calculate velocity vector

v_rot = (R1 * R2 * R3) * v_w';

% Define the cartesian coordinates

x = r_rot(1);

y = r_rot(2);

z = r_rot(3);

vx = v_rot(1);

vy = v_rot(2);

vz = v_rot(3);

end