代码：

function SunEarthMoon   % M函数文件load planets;   % 将planets.mat中的变量mass、position、velocity加载过来[sun, earth, moon] = deal(18, 3, 25);   % sun、earth、moon分别是18、3、25行
list = [sun, earth, moon];  % 1行3列矩阵
G = 6.67e-11; % gravitational constantdt = 24*3600;   % 作图的时间间隔为一天，每天有24*3600秒
N = length(list);   % N=3，三个天体
mass = mass(list);  % N行1列矩阵，N个天体的质量
position = position(list,:);    % N行3列矩阵，N个天体的坐标，坐标是1行3列的行向量，三个方向的分量
velocity = velocity(list,:);    % N行3列矩阵，N个天体的速度，速度是1行3列的行向量，三个方向的分量
h = plotplanets(position);  %作图函数for t = 1:365   % 图中总时间为一年，一年365天plotplanets(position,h);    % force = zeros(N,3); % N行3列零矩阵，一行表示某个天体在三个方向上的受力for i = 1 : N   % 遍历计算各天体间的万有引力。组合数C（3，2）Pi = position(i,:); % 某天体坐标Mi = mass(i);   % 某天体质量for j = (i+1):N     %the i+1 is to create diagonal Mj = mass(j);   % 另一天体质量Pj = position(j,:); % 另一天体坐标dr =  Pj - Pi;  % 两天体的相对，1行3列矩阵forceij = G*Mi*Mj./(norm(dr).^3).*dr;   % 两天体之间的力，1行3列的向量force(i,:) = force(i,:) + forceij;  % 规定正方向，将力计算进矩阵force(j,:) = force(j,:) - forceij;  % 反作用力与作用力方向相反，将力计算进矩阵% 上两行可替换为force([i,j],:) = force([i,j],:)+[forceij; -forceij];endendvelocity = velocity + force ./ repmat(mass,1,3)*dt;   % v=v+a*dt a=F/mposition = position + velocity*dt;  % r=r+v*dt
end % -------------------------------------------------------------------------function h = plotplanets(pos,h) % 
scale = 50;
total_planets = size(pos,1);
[sun, earth, moon] = deal(1, 2, 3);
radius = [50, 30, 20];
marker = {'.r', 'b.','m.'};
pos(moon,:) = pos(earth,:) + scale*(pos(moon,:)-pos(earth,:));
if nargin==1hold on; axis imageaxis( [-2 2 -2 2]*1e11 );for i = 1:total_planetsif any(i == [sun, earth, moon])h(i) = plot(pos(i,1),pos(i,2),marker{i},'markersize',radius(i));plot(pos(i,1), pos(i,2), marker{i}, 'markersize',5);elseh(i) = plot(pos(i,1), pos(i,2), 'k.', 'markersize', 20);plot(pos(i,1), pos(i,2), 'k.', 'markersize',5);endend
elsefor i = 1:total_planetsset(h(i), 'Xdata', pos(i,1)  , 'Ydata', pos(i,2)  )if any(i == [sun, earth, moon])plot(pos(i,1), pos(i,2), marker{i}, 'markersize',5);elseplot(pos(i,1), pos(i,2), 'k.', 'markersize',5);endenddrawnow
end

效果图