%
%   necessary input constants
%

m = 1.0;                      %  mass of the bar
g = 9.81;                   %  gravity
l = 1.0;                      %  length of the bar

t_final = 5;                %  time duration for calculation
                            %  initial time is always 0

%
%   Initial conditions for theta and theta_dot
%

theta0 = 1e-4;               %  initial displacement of the bar
                             %  initial disturbance in other words
theta_dot0 = 0;              %  initial angular velocity of the bar

%
%   Setup integration solver (ODE45)
%

ics = [theta0,theta_dot0];

%
%  Perform ODE45 solver
%

options = odeset('Refine',6,'RelTol',1e-7,'AbsTol',1e-10);
 
%  options 'Refine' was used to produce good-looking plots
%  default value of Refine was 4.
%  Refine doesn't apply if length(TSPAN) > 2.
%  scalar relative error tolerance 'RelTol' (1e-3 by default) and
%  vector of absolute error tolerances 'AbsTol' (all components 1e-6 by
%  default).

[t,x] = ode45('pendulum_eom',[0 t_final],ics,options);

Theta = x(:,1);
Theta_dot = x(:,2);
subplot(2,1,1);
plot(t,Theta,'g');title('Angular Displacement vs. Time');
ylabel('\theta (rad)');xlabel('time (s)');
subplot(2,1,2);
plot(t,Theta_dot,'r');title('Angular Velocity vs. Time');
ylabel('d\theta/dt (rad/s)');xlabel('time (s)');

Pendulum problem (subprogram source must be saved as pendulum_eom.m)

function xdot = pendulum_eom(t,x)

%
%   Constants and input variables
%

g = 9.81;            %  gravity
l = 1;               %  length of the bar

x_dot1 = x(2);
x_dot2 = -(3*g)/(2*l)*sin(x(1));

xdot = [x_dot1; x_dot2];