% ch2example6main.m
clear;
v0=0; y0=1;     % 球的初始狀態(tài)
x_state=[v0,y0];% 將初始狀態(tài)賦值到狀態(tài)變量中
dt=0.01;        % 仿真步進
t=0:dt:5;       % 仿真時間序列
K=0.85;         % 碰撞衰減系數(shù)
for k=1:length(t)   % 仿真開始，每次循環(huán)向前推進一個仿真步進dt
    x(k,:)=x_state; % 記錄并保存當(dāng)前狀態(tài)的計算結(jié)果
    [t_out, x_out]=ode45('ch2example6statefun',[t(k),t(k)+dt],x_state);
    % 計算下一個時刻的新狀態(tài)
    % 可換用ode23、ode113、ode23t、ode15s、ode23s以及ode23tb求解器
    x_state=x_out(length(x_out),:);     % 更新狀態(tài)
    if (x_state(2)<=0) & (x_state(1)<0) % 當(dāng)速度為負(fù)(球向下運動)且已經(jīng)接觸碰撞面
        x_state(1)=-K*x_state(1);       % 處理碰撞瞬間情況：速度反向并衰減K
    end
    % 動畫作圖：顯示小球彈跳過程
    y=x_state(2); % 小球當(dāng)前位置
    subplot(2,1,1);plot(0,y,'o'); 		 
    axis([-2 2 -0.1 1]); 	        % 坐標(biāo)范圍固定	 
    set(gcf,'DoubleBuffer','on'); 	% 雙緩沖避免作圖閃爍	 
    drawnow; % 立即顯示作圖
end
% 仿真結(jié)束。最后輸出計算時間序列上的結(jié)果：隨時間變化的速度和位移曲線
subplot(2,1,2);plotyy(t,x(:,1),t,x(:,2),'plot');