function pendulum2
% pendulum with damping constant c>=0  
% theta(t) = y(t) is angle from bottom position at time t
% ODE is  y'' + sin(y) = -c*y'

[t,y] = ode23(@f,0:.01:80,[pi-.1;0]);   % want output for t=0, .01, .02, ... , 80

theta = y(:,1);                         % angle theta for t=0, .01, .02, ... , 80
thetadot = y(:,2);                      % derivative of theta(t)
plot(sin(theta),-cos(theta)); 
axis([-1.1 1.1 -1.1 1.1]); axis equal; hold on
xlabel('x'); ylabel('y')
title('trajectory of bob of pendulum')

disp('hit any key for animation'); pause
comet(sin(theta),-cos(theta));   hold off

figure
plot(theta,thetadot); axis equal; hold on
xlabel('\theta'); ylabel('\theta''')
title('phase plane');

disp('hit any key for animation'); pause
comet(theta,thetadot);  hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function yp = f(t,y)  
c = 0.1;                % with damping
yp = [ y(2) ; -sin(y(1)) - c*y(2) ];
