function dydt=damped(t,y)
% DAMPED(T,Y) where T is scalar and Y=[Y1;Y2;Y3;Y4].
% Sample file suitable for (numerically) integrating a
% nonlinear oscillator and its linearization.
% Y1 is the position, Y2 the velocity of the nonlinear
% oscillator.
% Y3 is the position and Y4 is the velocity of the linear
% oscillator.
% The "damping constants" p1 and p3, as well as the 
% "inertia" q1 and q3 are set inside this file.
%
% Note that upon the substitutions y1=y and y2=y',
% and y3=y and y4=y', respectively, the 2nd order DE
% are written as the systems of 1st order DEs
% y''+p2*y'+q2*sin(y)=0 becomes y1'=y2, y2'=-p1*y2-q1*sin(y3)
% y''+p2*y'+q2*y=0 becomes y3'=y4, y4'=-p2*y4-q2*y3
%
% Original version April 1999.   UPDATE May 1999: 
% The signs of the p's and q's have been changed!
%
% All rights reserved. Matthias Kawski. May 1999.

p1=0.2;
p3=p1;
q1=3;
q3=q1;
dydt=[y(2);... 
      -q1*sin(y(1))-p1*y(2);...
      y(4);...
      -q3*y(3)-p1*y(4)];