% Constants
m=15;
l=2;
g=9.81;
d=l/2;
Io=(1/3)*m*(l^2);
wn=sqrt((m*g*d)/Io);
% function handle for linearized ode
Vib_lin = @(~,theta,wn) [theta(2); -(wn^2).*(theta(1))];
% function handle for nonlinearized ode
Vib_nonlin = @(~,theta,wn) [theta(2); -(wn^2).*sin(theta(1))];
% input data for ode45
theta0=5;
thetadot0=0;
t=[0 10];
% solve linearized part
[t_lin,theta_lin]=ode45(@(~,theta)Vib_lin(t,theta,wn),t,[theta0 thetadot0]);
% solve nonlinearized part
[t_nonlin,theta_nonlin]=ode45(@(~,theta)Vib_nonlin(t,theta,wn),t,[theta0 thetadot0]);
% plot results
figure (1);
plot(t_lin,theta_lin(:,1),t_nonlin,theta_nonlin(:,1));