% 複数の固定ステップ幅で微分方程式を解く

global angularfrequency;
angularfrequency = 2*pi;

timefinal = 5.00;
qinit = [2.00;0.00];

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

for timestep = [0.1 0.01 0.001]
    interval = 0.00:timestep:timefinal;
    [time,q]=ode45(@simple_vibration,interval,qinit);
    plot(time,q(:,1),'-');
    ylim([-2.5,2.5]);
    xlabel('time');
    ylabel('angle');
    grid on;
    filename = strcat('simple-vibration-', num2str(timestep), '.eps');
    saveas(gcf, filename, 'eps');
    fprintf("固定ステップ幅 %f\n", timestep);
    fprintf("一時停止：何かのキーを押してください\n");
    pause;
end

% 可変ステップ幅で微分方程式を解く

interval = [0.00, timefinal];
[time,q]=ode45(@simple_vibration,interval,qinit);
plot(time,q(:,1),'-');
ylim([-2.5,2.5]);
xlabel('time');
ylabel('angle');
grid on;
saveas(gcf, 'simple-vibration-adaptive.eps', 'eps');
fprintf("可変ステップ幅\n");
fprintf("一時停止：何かのキーを押してください\n");
pause;

s = size(time);
n = s(1);
dt = time(2:n) - time(1:n-1);
plot(time(1:n-1),dt);
ylim([0,0.042]);
xlabel('time');
ylabel('time step');
yticks([0 0.010 0.020 0.030 0.040]);
grid on;
saveas(gcf, 'simple-vibration-interval.eps', 'eps');
fprintf("各時刻におけるステップ幅\n");
fprintf("一時停止：何かのキーを押してください\n");
pause;

energy = (1/2)*q(:,2).^2 + (angularfrequency )^2 * (1/2)*q(:,1).^2;
plot(time,energy,'-');
grid on;
fprintf("各時刻におけるエネルギー\n");
fprintf("一時停止：何かのキーを押してください\n");
pause;
