ode45 to solve second order ODE

Hey guys,

I'm trying to use Matlab to solve a second-order ODE that's basically an equation of motion derived from the Lagrangian equation. I'm supposed to solve it using various initial conditions, then theta(0)=10 degrees, theta(0)=20 degrees, theta(0)=30 degrees, theta(0)=40 degrees, all whereby thetadot(0) (basically diff(theta,t)) = 0.

Below is my code:

syms m g p theta(t) r thetadot I Theta
%I = m*r^2*(3/2*8/(3*pi));
omega=diff(theta,t)
h = p*(1-cos(theta));
T=0.5*I*(omega^2);
V=m*g*h;
%d2y = diff(theta,t,2)
L = T-V
L1 = subs(L,omega,thetadot);
L1diff = diff(L1,thetadot);
L2 = subs(L,theta,Theta)
L2new=diff(L2,Theta)
L3 = subs(L1diff,thetadot,omega)
L3diff = diff(L3,t)
EOM = L3-L2 == 0
cond2 = 0; %dy(0)
cond1 = 0.174; %theta(0)

I have been able to get my equation of motion and I was given such conditions, but I have no idea how to solve a second-order ODE with ode45, and the documentation is rather abstract. Any help would be greatly appreciated. Thanks all!

[Star Strider]

I can’t figure out what your differential equation is, even after running your code.

Use the odeToVectorField function to create a vector field from your second-order ODE, then matlabFunction to convert the first output of odeToVectorField to an anonymous function that you can use with ode45.

[Isaac Tee]

My apologies. I've made a few errors in my ODE. Here is my latest code:

syms m g p theta(t) r thetadot I Theta x(t) y(t)
%I = m*r^2*(3/2*8/(3*pi));
omega=diff(theta,t)
x = p*sin(theta);
h = p*(1-cos(theta));
y = h;
xdot = diff(x,t)
ydot = diff(y,t)
T=0.5*m*(xdot^2+ydot^2)+0.5*I*(omega^2)
V=m*g*h
%d2y = diff(theta,t,2)
L = T-V
L1 = subs(L,omega,thetadot);
L1diff = diff(L1,thetadot)
L2 = subs(L,theta,Theta)
L2new=diff(L2,Theta);
L2final = subs(L2new,Theta,theta)
L3 = subs(L1diff,thetadot,omega);
L3diff = diff(L3,t)
EOM = L3diff-L2final == 0
% The ODE should be (m*(2*p^2*cos(theta(t))^2*diff(theta(t), t, t) + 2*p^2*sin(theta(t))^2*diff(theta(t), t, t)))/2 + I*diff(theta(t), t, t) + g*m*p*sin(theta(t)) == 0
cond2 = omega(0) == 0.174; %dthetadt(0)
cond1 = theta(0)==0;
conds = [cond1,cond2];