syms y(t) theta(t)
g = 9.8;
m = 1;
k = 100;
Ln = .5;
y0 = .2;
ydot0 = 0;
theta0 = 1;
thetadot0 = 0;
initialVals = [.2,1,0,0];
eqn1 = diff(y,t,2) == g*cosd(theta) + (Ln+y)*diff(theta,t)^2 - k*y/m;
eqn2 = diff(theta,t,2) == (-g*sind(theta)-2*diff(theta,t)*diff(y,t))/(Ln+y);
vars = [y(t);theta(t)];
eqnSystem = odeToVectorField([eqn1,eqn2]);
eqnSystem = matlabFunction(eqnSystem,'vars', {'t','Y'});
thetaVals = ode45(eqnSystem,[0,16],initialVals);
tVec = linspace(0,16);
thetaVec = deval(thetaVals,tVec,1);
plot(tVec,thetaVec);
title('When Ln = .5m')
xlabel('Time (seconds')
ylabel('Pendulum angle (degrees)')