% 硬い常微分方程式を解く
% 衝突の運動方程式を数値的に解く

global mass;
global spring;
global grav;

mass = 1.0;
spring = 100.0;
grav = 9.8;

interval = [0, 50]; % 可変ステップ
qinit = [100; 0]; % 初期値
options = odeset('AbsTol',1e-12, 'MaxStep',0.1);
%[time, q] = ode45(@collision, interval, qinit, options);
[time, q] = ode23tb(@collision, interval, qinit, options);

set(0,'defaultAxesFontSize',16);
set(0,'defaultTextFontSize',16);

plot(time,q(:,1));
xlabel('time');
ylabel('position');
%legend('x');
grid on;
saveas(gcf, 'collision-position.eps', 'eps');
fprintf("t-x グラフ\n");
fprintf("一時停止：何かのキーを押してください\n");
pause;

plot(time,q(:,2));
xlabel('time');
ylabel('velocity');
%legend('v');
grid on;
saveas(gcf, 'collision-velocity.eps', 'eps');
fprintf("t-v グラフ\n");
fprintf("一時停止：何かのキーを押してください\n");
pause;

s = size(time);
n = s(1);
dt = time(2:n) - time(1:n-1);
plot(time(1:n-1),dt);
xlabel('time');
ylabel('time step');
grid on;
saveas(gcf, 'collision-interval.eps', 'eps');
fprintf("t-ステップ幅\n");
fprintf("一時停止：何かのキーを押してください\n");
pause;
