function dy=Rayleigh(t,y)
% Rayleigh oscillator

u1=1;
u2=878e-6;
a1=1.542;
a2=0.5;
b=10e-6;
c=2.377e-6;
omega=60.3e3;
dy=zeros(2,1);

% Nonlinearity
T=-a1*(-y(1))+a2*(-y(1))^3;
% System
dy(1)=(1/c)*(y(2)+T);
dy(2)=(1/b)*(-y(1)-u1-u2*sin(omega*t));


% "Complex and Chaotic Nonlinear Dynamics. 
% Advances in Economics and Finance, 
% Mathematics and Statistics" 
% T.Vialar, Springer 2009
% Copyright(c).