function [t,phi] = pendulum_RK4(T,dt,phi0)

% function [t,phi] = pendulum_RK4(T,dt,phi0)
%
% 4th order Runge-Kutta applied to nonlinear pendulum eqn:
%   dphi(1)/dt = phi(2)                 
%   dphi(2)/dt = -sin(phi(1))
% Arguments:
%   T     final integration time 
%   dt    timestep  
%   phi0  initial phi at time t=0
% Outputs:
%   t     Vector of times 0:dt:T of length nt
%   phi   Matrix phi(nt,2) of solutions at the times t. 
  
F = inline('[phi(2) -sin(phi(1))]','phi');
nt = round(T/dt);
phi = zeros(nt,2);
phi(1,:) = phi0;
for n = 1:nt
  q1 = dt*F(phi(n,:));
  q2 = dt*F(phi(n,:)+0.5*q1);
  q3 = dt*F(phi(n,:)+0.5*q2);
  q4 = dt*F(phi(n,:)+q3);
  phi(n+1,:) = phi(n,:) + (q1+2*q2+2*q3+q4)/6;
end

t = (0:nt)*dt;