% pendulum_test.m
% Generating plots for non-linear pendulum

m = 1; % mass in kg
L= 1;  % length in meters
g = 9.81; % acceleration due to gravity

tspan = [0,20]; % range of times from 0 to 20 seconds

c = input('Damping coefficient: ');
theta = input('Initial angle in degrees: ');
vel = input('Angular velocity in degrees/sec: '); 

init = (pi/180)*[theta,vel];

options = odeset('RelTol' ,1e-6, 'AbsTol', 1e-6);

[t,x] = ode45(@(t,x) pendulum(t,x,m,L,g,c), ...
        tspan, init, options);

% phase plane: angle on horizontal axis
% angular velocity on vertical axis

subplot(2,2,1);
plot(x(:,1), x(:,2), '-');

% Angle as a function of time

subplot(2,2,2);
plot(t, x(:,1), '-');

% Angular velocity as a function of time

subplot(2,2,3);
plot(t, x(:,2), '-');
shg;