%  Script to plot RK4 method with dt = 4*pi/10 and pi/10 applied to
%  the nonlinear pendulum equation 
%    dphi(1)/dt = phi(2)
%    dphi(2)/dt = -sin(phi(1))
%  with IC of pendulum starting at rest 1 radian from the vertical:
%    phi(1) = 1, phi(2) = 0
%  Plots show phi(1) over roughly the 1-2 and the 49-50 pendulum
%  periods for each dt.

T = 100*pi;
phi0 = [1 0];
trange1 = [0 2]*2*pi;
trange2 = [48 50]*2*pi;

dt = pi/10;
[t,phi] = pendulum_RK4(T,dt,phi0);
subplot(2,2,1)
irange1 = 1 + round(trange1(1)/dt):round(trange1(2)/dt);
plot(t(irange1)/(2*pi), phi(irange1,1))
title('RK4, \Delta t = \pi/10')
xlabel('t/2\pi')
ylabel('\phi_1')
subplot(2,2,2)
irange2 = 1 + round(trange2(1)/dt):round(trange2(2)/dt);
plot(t(irange2)/(2*pi), phi(irange2,1))
xlabel('t/2\pi')
ylabel('\phi_1')


dt = 4*pi/10;
[t,phi] = pendulum_RK4(T,dt,phi0);
subplot(2,2,3)
irange1 = 1 + round(trange1(1)/dt):round(trange1(2)/dt);
plot(t(irange1)/(2*pi), phi(irange1,1))
title('RK4, \Delta t = 4\pi/10')
xlabel('t/2\pi')
ylabel('\phi_1')
subplot(2,2,4)
irange2 = 1 + round(trange2(1)/dt):round(trange2(2)/dt);
plot(t(irange2)/(2*pi), phi(irange2,1))
xlabel('t/2\pi')
ylabel('\phi_1')
