%% integrate the differential equation of a pendulum

function pendulum()
  clear;
  a=4; % length of first arm
  g=9.81; % acceleration of gravity
  t=0:.05:10; % times to display
  q=[2 0]; % initial conditions

  [t,q]=ode45( @fn, t, q );
  %%[t,q]=ode45( @(t,q)[q(2);-g/a*sin(q(1))], t, q );
  %%plot(t,q(:,1))

  %% animation of motion
  clf;
  axis([-5 5 -5 5],'manual','equal');
  arm=line([0,0],[0,-a],'Marker','o');
  l=line([5,-5,-5,5,5],[5,5,-5,-5,5]);
  for i=1:size(t)
    set(arm, 'XData', [0, a*sin(q(i,1))]);
    set(arm, 'YData', [0,-a*cos(q(i,1))]);
    drawnow; pause(.05);
  end
end

function d=fn(t,q)
  a=4; % length of first arm
  g=9.81; % acceleration of gravity
  d=[q(2);-g/a*sin(q(1))];
end
