ODE

function [dtheta_dt] = ode_func(t,theta,b,g,l,m)

theta1=theta(1)
theta2=theta(2)
dtheta1_dt=theta2;
dtheta2_dt=-(b/m)*theta2-(g/l)*sin(theta1);
dtheta_dt=[dtheta1_dt;dtheta2_dt];
end

clear all
close all
clc

b=0.5;
g=9.81;
l=1;
m=1;
% initial condition
theta_0=[0;3];
% time points
tspan=linspace(0,20,500);
% solve ODE
[t,results]=ode45(@(t,theta) ode_func(t,theta,b,g,l,m),tspan,theta_0);
plot(t,results(:,1))
hold on
plot(t,results(:,2))
xlabel(\'time\')
ylabel(\'plot\')