t = 10;
L = 10;
g = 9.81;
m = 1;
syms m a g theta(t)
eqn = m*a == -m*g*sin(theta);
syms L
eqn = subs(eqn,a,L*diff(theta,2));
eqn = isolate(eqn,diff(theta,2));
syms omega_0
eqn = subs(eqn,g/L,omega_0^2);
syms x
approx = taylor(sin(x),x,'Order',2);
approx = subs(approx,x,theta(t));
eqnLinear = subs(eqn,sin(theta(t)),approx);
syms theta_0 theta_t0
theta_t = diff(theta);
cond = [theta(0) == theta_0, theta_t(0) == theta_t0];
assume(omega_0,'real')
thetaSol(t) = dsolve(eqnLinear,cond)
omega_0Value = sqrt(g/L);
T = 2*pi/omega_0Value;
theta_0Value = 0.1*pi;
theta_t0Value = 0;
vars = [omega_0 theta_0 theta_t0];
values = [omega_0Value theta_0Value theta_t0Value];
thetaSolPlot = subs(thetaSol,vars,values);