高手救命啊！

newton_grc

function hit
clear all;
clc;
y0 = [1;0;0;1;
1;0;1;-1;] %粒子的初位置和初速度


tn=[0];
for i=1:2%循环两次

opts = odeset('events',@events1);%控制方解程
[t,y,te]=ode45(@f,[tn(1) Inf],y0,opts);
comet(y(:,3),y(:,1));
plot(y(:,3),y(:,1)); %画出粒子1的轨迹
hold on;
comet(y(:,7),y(:,5));
plot(y(:,7),y(:,5),'m'); %画出粒子2的轨迹
hold on;

if i<3
tn=[te(i)]
end
if i>=3
tn=[te(2)]
end %每次碰撞结束后从新设定方程的时间起点

y01=[]
for n=1:4:8
if y(end,n)<0.0001
y01=[y01;y(end,n);-y(end,n+1)]%当y=0时即碰到下底时y方向速度改变
else
y01=[y01;y(end,n);y(end,n+1)]
end
if abs(y(end,n+2)-2)<0.0001|abs(y(end,n+2)+2)<0.0001
y01=[y01;y(end,n+2);-y(end,n+3)]%当x=-2时即碰到两个侧壁时x方向速度改变
else
y01=[y01;y(end,n+2);y(end,n+3)]
end
for m=(n+4):4:8

if (abs(y(end,n)-y(end,m))+abs(y(end,n+2)-y(end,m+2)))<0.001%当两个粒子相碰的时候

linshiy=y01(n+1);
y01(n+1)=y01(m);
y01(n+1)=linshiy;
linshix=y01(n+3);
y01(n+3)=y01(m+3);
y01(m+3)=linshix;
%交换速度
end
end
end
y0=y01
end
xlabel('x'); ylabel('y');title('重力场中两粒子相碰')

function ydot=f(t,y)
odes=[];
for i=2:4:8
odes=[odes;y(i);-1;y(i+2);0];
end
ydot=odes;


function [value,isterminal,direction] = events1(t,y)
valuem=[]
for i=1:4:8
valuem=[valuem,y(i),y(i+2)-2,y(i+2)+2]
for n=i+4:4:8
valuem=[valuem,abs(y(i)-y(n))+abs(y(i+2)-y(n+2))]
end
end


value=valuem
isterminal =ones(1,7);
direction = zeros(1,7);

程序中红色部分表示连个例子相碰的判断条件以及相碰后的速度及位移。
存在两个问题：
一、由于ode45计算步长的缘故两个粒子相碰时并没有停止计算
二、方程在有些初始数值时会出现错误。比如现在给的这组值，急需高手救命！

[ 本帖最后由 newton_grc 于 2006-7-22 20:46 编辑 ]