function main()
global m mu g

m=1; % 小球质量
mu=0.3; % 摩擦系数
V0=100; % 初始速度
theta=pi/6; % 初始角度
g=9.8; % 重力加速
k=0.8; % 撞击系数
H=10; % 初始高度

tspan=[0 10]; % 求解时间范围
x0=[0,V0*cos(theta),H,V0*sin(theta)]; % 微分方程初值

[t,x]=ode45(@odefun,tspan,x0);
X=x(:,1);
Vx=x(:,2);
Y=x(:,3);
Vy=x(:,4);
figure
plot(t,x)
legend('X','Vx','Y','Vy')
xlabel('时间')
ylabel('变量')
grid on
figure
plot(X,Y)
xlabel('水平方向')
ylabel('竖直方向')
title('小球运动轨迹')
grid on

function dx=odefun(t,x)
% 微分方程描述函数

global m mu g
% X=x(1);
Vx=x(2);
% Y=x(3);
Vy=x(4);

% 运动微分方程
dx=[Vx
    -mu*Vx/m       % dVx/dt
    Vy
-mu*Vy/m-g];   % dVy/dt

function main
global m mu g

m=1; % 小球质量
mu=0.3; % 摩擦系数
V0=100; % 初始速度
theta=pi/6; % 初始角度
g=9.8; % 重力加速
k=0.8; % 撞击系数
H=10; % 初始高度

tspan=[0 10]; % 求解时间范围
x0=[0,V0*cos(theta),H,V0*sin(theta)]; % 微分方程初值

options=odeset('events',@events); % 设置小球过零点（撞地）检测条件
[t,x,te,xe,ie]=ode45(@odefun,tspan,x0,options);

X=x(:,1);
Vx=x(:,2);
Y=x(:,3);
Vy=x(:,4);
figure
plot(t,x)
legend('X','Vx','Y','Vy')
xlabel('时间')
ylabel('变量')
title('参数随时间变化图形')
grid on
figure
plot(X,Y)
xlabel('水平方向')
ylabel('竖直方向')
title('小球运动轨迹图形')
grid on

function [value,isterminal,direction]=events(t,x)
% 事件检查函数，此时需要做的是过零点检测
Y=x(3);

% ode45函数自动检查当value=0是否成立
% 在此问题上我们要求检测Y=0的点，因此设置value=Y
% 如果我们要检测Y=2，那么就设置value=Y-2
value=Y;
% 检测到指定条件时，是否终止ode45函数的运行
% 1表示终止，0表示继续
% 在我们这个问题上，我们只要检测到零点时就停止程序
isterminal=1;
% value过零点检测的方向
% -1表示由正到负，+1表示由负到正
% 对于我们这个问题，当然是由正到负
direction=-1;

function main
global m mu g

m=1; % 小球质量
mu=0.25; % 摩擦系数
V0=100; % 初始速度
theta=pi/6; % 初始角度
g=9.8; % 重力加速
k=0.9; % 撞击系数
H=10; % 初始高度

tspan=[0 100]; % 求解时间范围
x0=[0,V0*cos(theta),H,V0*sin(theta)]; % 微分方程初值
options=odeset('events',@events); % 设置小球过零点（撞地）检测条件

n=10; % 反弹次数
% 初始化必要的参数
tout=tspan(1); % 记录完整过程的时间
xout=x0; % 记录完整过程的参数
teout=[]; % 记录撞击点时刻
xeout=[]; % 记录撞击点参数

for i=1:n
    [t,x,te,xe]=ode45(@odefun,tspan,x0,options);
    % 保存该阶段的参数
    tout=[tout;t(2:end)];
    xout=[xout;x(2:end,:)];
    teout=[teout;te];
    xeout=[xeout;xe];
   
    % 修改方程求解区间
    tspan(1)=te;
    % 修改初值方程
    x0(1)=xe(1); % 水平位移
    x0(2)=xe(2)*k; % 水平速度，注意乘以撞击系数
    x0(3)=0; % 竖直位移
    x0(4)=-xe(4)*k; % 竖直速度，注意乘以撞击系数，并且反弹
end

X=xout(:,1);
Vx=xout(:,2);
Y=xout(:,3);
Vy=xout(:,4);
figure
plot(tout,xout)
hold on
plot(teout,xeout,'ro','MarkerFaceColor','red')
legend('X','Vx','Y','Vy')
xlabel('时间')
ylabel('变量')
title('参数随时间变化图形')
grid on
figure
plot(X,Y)
hold on
plot(xeout(:,1),xeout(:,3),'ro','MarkerFaceColor','red')
xlabel('水平方向')
ylabel('竖直方向')
title('小球运动轨迹图形')
grid on