function vprime = odeplay(t,vin) global r; global a; global b; global k; x = vin(1); y = vin(2); z = vin(3); % compute derivatives for damped pendulum, % from K Tung p. 227 xprime = y; yprime = -b*sin(x) - k*y + a*cos(z); zprime = r; % map back to vector format vprime=[xprime;yprime;zprime]; end