function ydot = pendrhs(t,y,par)
%PENDRHS Compute right-hand side of double pendulum model
%   Inputs are time (ignored), a state vector consisting of two angles and
%   two generalized momenta, and a vector of parameters.

theta1 = y(1);
theta2 = y(2);
p1 = y(3);
p2 = y(4);
ydot = zeros(4,1);

temp = par(2)*cos(theta1-theta2);
denom = par(3)*par(4) - temp^2;

ydot(1) = (par(4)*p1 - temp*p2)/denom;
ydot(2) = (par(3)*p2 - temp*p1)/denom;
ydot(3) = -par(5)*par(1)*sin(theta1) - par(2)*sin(theta1-theta2)*ydot(1)*ydot(2);
ydot(4) = -par(5)*par(2)*sin(theta2) + par(2)*sin(theta1-theta2)*ydot(1)*ydot(2);

end

