matlab解微分方程组除了初值全是NaN

代码如下：
main.m:

clear all;
close all;
%ode45求解
tspan=[0:50:50000];%求解区间
a0=[0,0,0,0,0,0];%初值[x,x',y,y',z,z']
[t,a]=ode45('odefun',tspan,a0);%使用ode45求解
plot3(a(:,1),a(:,3),a(:,5),'*');

odefun.m:

function da = odefun(t,a)
% 定义
% a(1)=x
% a(2)=x'
% 定义
% a(3)=y
% a(4)=y'
% 定义
% a(5)=z
% a(6)=z'
rp=[476333.10266411;0.00658345;4478.74238999];%探测器相对于地球质心的位置矢量
rs=[10;120;1303];%太阳相对于地球质心的位置矢量
rm=[1;2;3];%月球相对于地球质心的位置矢量

rpm=norm(rm-rp);
rps=norm(rs-rp);
rms=norm(rs-rm);

us= 1.3272e+20;%太阳的引力常数
ue= 3.9802e+14;%地球的引力常数
um= 4.8990e+12;%月球的引力常数

hm=cross(rm,rm');%月球轨道的瞬时角动量

%地心会合坐标系相对于惯性坐标系的角速度
wx=norm(rm)/norm(hm)^2*(us/rms^3-us/norm(rs))*dot(rs,hm);
wy=0;
wz=norm(hm)/norm(rm)^2;
w=[wx;wy;wz];

rsx=dot(rs,rm/norm(rm));
rsy=dot(rs,cross(hm,rm)/(norm(rm)*norm(hm)));
rsz=dot(rs,hm/norm(hm));

da=zeros(6,1);%初始化
da(1)=a(2);
da(2)=2*wz*a(4)+wz^2*a(1)-wx*wz*a(5)+wz'*a(3)-(ue/norm(rp)^3+um/norm(rpm)^3+us/norm(rps)^3)*a(1)+(um/norm(rpm)^3-um/norm(rm)^3)*norm(rm)+(us/norm(rps)^3-us/norm(rs)^3)*norm(rsx);
da(3)=a(4);
da(4)=2*wx*a(6)-2*wz*a(2)-(wx^2+wz^2)*a(3)+wx'*a(5)-wz'*a(1)-(ue/norm(rp)^3+um/norm(rpm)^3+us/norm(rps)^3)*a(3)+(us/norm(rps)^3-us/norm(rs)^3)*norm(rsy);
da(5)=a(6);
da(6)=-2*wx*a(4)+wx^2*a(5)-wx*wz*a(1)+wx'*a(3)-(ue/norm(rp)^3+um/norm(rpm)^3+us/norm(rps)^3)*a(5)+(us/norm(rps)^3-us/norm(rs)^3)*norm(rsz);
end

结果如图：

你好

hm=cross(rm,rm');%月球轨道的瞬时角动量

这个语句有问题。两个相同矢量的叉乘结果是0矢量，任何数除以0都是NaN，建议仔细看看