%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ダフィングの方程式
% d^2x/dt^2 + 2γ dx/dt + ω0^2 x + ε x^3 = A cos(ω t)
% の定義
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dxdt = Duffing(t, x, gamma, omega0, epsilon, A, omega)
%%%%%% 引数 %%%%%%
% t 時間
% x 状態量
%    x(1)がx(振動変位)
%    x(2)がdx/dt(振動速度)を表している
%    を表している
% gamma γ，2*γ = b/m(制動係数)
% omega0 ω0 = sqrt(k/m)(固有角振動数)
% epsilon ε = β/m
% A 加振力の振幅
% omega ω(加振力の角振動数)
%%%%%% 戻り値 %%%%%%
% dxdt xの1階微分値
%    dxdt(1)がdx/dt(= x(2))
%    dxdt(2)がd^2x/dt^2(= - 2γ dx/dt - ω0^2 x - ε x^3 + A cos(ω t))
%    を表している



dxdt = [x(2); -2*gamma*x(2) - omega0^2*x(1) - epsilon*x(1)^3 + A*cos(omega*t)];
