请用matlab语言写一份代码来绘制三维哈雷彗星运行轨迹，要求代码适应R2021a版本的matlab，

好的，这里是一份用MATLAB语言编写的绘制哈雷彗星运行轨迹的代码，适用于R2021a版本的MATLAB：

% 定义哈雷彗星的轨道参数
a = 17.834e10; % 半长轴（m）
e = 0.96714; % 离心率
i = 162.262; % 轨道倾角（°）
Omega = 58.42; % 升交点赤经（°）
omega = 111.33; % 近地点幅角（°）
M0 = 38.38; % 平均近点角（°）
T = 76.1*365*24*3600; % 周期（s）

% 定义计算时间范围
t_start = 0;
t_end = 2*T;
dt = 3600;

% 计算哈雷彗星在每个时刻的位置
t = t_start:dt:t_end;
n = 2*pi/T;
M = M0 + n*t;
E = M;
for i = 1:length(t)
    while abs(E(i) - e*sin(E(i)) - M(i)) > 1e-8
        E(i) = E(i) - (E(i) - e*sin(E(i)) - M(i))/(1 - e*cos(E(i)));
    end
end
nu = 2*atan(sqrt((1+e)/(1-e)).*tan(E/2));
r = a*(1-e^2)./(1+e*cos(nu));
x = r.*cos(nu);
y = r.*sin(nu);
z = zeros(size(x));
[xx,yy,zz] = sph2cart(deg2rad(Omega),deg2rad(i),deg2rad(omega));
R = [cosd(Omega)*cosd(omega)-sind(Omega)*cosd(i)*sind(omega) -cosd(Omega)*sind(omega)-sind(Omega)*cosd(i)*cosd(omega) sind(Omega)*sind(i); ...
     sind(Omega)*cosd(omega)+cosd(Omega)*cosd(i)*sind(omega) -sind(Omega)*sind(omega)+cosd(Omega)*cosd(i)*cosd(omega) -cosd(Omega)*sind(i); ...
     sind(i)*sind(omega) cosd(i)*sind(omega) cosd(i)*cosd(omega)];
pos = R*[x; y; z];
x = pos(1,:);
y = pos(2,:);
z = pos(3,:);

% 绘制哈雷彗星运行轨迹
figure;
plot3(x, y, z);
title('Halley Comet Orbit');
xlabel('x (m)');
ylabel('y (m)');
zlabel('z (m)');
axis equal;

此代码将计算哈雷彗星在给定时间范围内的位置，并绘制其三维运行轨迹。您可以在MATLAB中运行此代码，生成哈雷彗星运行轨迹的图形。