% ch2example5prg1.m
dt=0.01;            % 仿真步進
T=15;               % 仿真時間長度
t=0:dt:T;           % 仿真計算時間序列
g=9.8;              % 重力加速度
L=1;                % 擺線長度
m=10;               % 擺錘質量
k=5;                % 空氣阻力比例系數
theta0=3.1;         % 初始擺角設置
v0=0;               % 初始擺速設置
x0=[v0;theta0];     % 初始狀態賦值
par=[g;k;m;L];      % 系統參數賦值
% 以歐拉算法計算并作圖
[t_out, x_out]=eulerode('pendulumstateeq',t', x0, [], par);
plot(t_out, x_out(:,1),'-.k');hold on;
% 以ode45算法計算并作圖對比
[t_out, x_out]=ode45('pendulumstateeq',t', x0, [], par);
plot(t_out, x_out(:,1),'-k');
xlabel('時間 t ');ylabel('線速度 m/s');
legend('歐拉算法','ode45算法');