function [y,dy] = dpend(t,y)
%Pendulum differential equations
% g [m/s^2] acceleration due to gravity
% length [m] pendulum length
% Izzi Urieli - Jan 21, 2002 (updated 12/08/07)
global g length
THETA = 1; % index of angle [radians]
OMEGA = 2; % index of angular velocity [rads/sec]
HEIGHT = 3; % index of pendulum mass height [m]
dy(THETA) = y(OMEGA);
dy(OMEGA) = -(g/length)*sin(y(THETA));
y(HEIGHT) = length*(1 - cos(y(THETA)));