Trying to solve systems of ODEs and graph answer

Trying to solve system of ODEs using other working code but getting error on line with "matlabFunction" I beleave the error has somethig to do with 'vars' but I can not figure out what. ODEs equation of motrions are correct and code works, error is just with graphing.

Original Code

clear;clc;

syms theta_1(t) theta_2(t) L_1 L_2 m_1 m_2 g

x_1 = L_1*sin(theta_1);
y_1 = -L_1*cos(theta_1);
x_2 = x_1 + L_2*sin(theta_2);
y_2 = y_1 - L_2*cos(theta_2);

vx_1 = diff(x_1);
vy_1 = diff(y_1);
vx_2 = diff(x_2);
vy_2 = diff(y_2);

ax_1 = diff(vx_1);
ay_1 = diff(vy_1);
ax_2 = diff(vx_2);
ay_2 = diff(vy_2);

syms T_1 T_2

eqx_1 = m_1*ax_1(t) == -T_1*sin(theta_1(t)) + T_2*sin(theta_2(t))
eqy_1 = m_1*ay_1(t) == T_1*cos(theta_1(t)) - T_2*cos(theta_2(t)) - m_1*g
eqx_2 = m_2*ax_2(t) == -T_2*sin(theta_2(t))
eqy_2 = m_2*ay_2(t) == T_2*cos(theta_2(t)) - m_2*g

Tension = solve([eqx_1 eqy_1],[T_1 T_2]);

eqRed_1 = subs(eqx_2,[T_1 T_2],[Tension.T_1 Tension.T_2]);
eqRed_2 = subs(eqy_2,[T_1 T_2],[Tension.T_1 Tension.T_2]);

L_1 = 1;
L_2 = 1.5;
m_1 = 2;
m_2 = 1;
g = 9.8;
eqn_1 = subs(eqRed_1)

eqn_2 = subs(eqRed_2)

[V,S] = odeToVectorField(eqn_1,eqn_2)

S

M = matlabFunction(V,'vars',{'t','Y'})

initCond = [pi/4 0 pi/6 0];
sols = ode45(M,[0 10],initCond);

plot(sols.x,sols.y)
legend('\theta_2','d\theta_2/dt','\theta_1','d\theta_1/dt')
title('Solutions of State Variables')
xlabel('Time (s)')
ylabel('Solutions (rad or rad/s)')

New codewith matlabFunction error

clear; clc;

% Initalize Variables
syms theta1(t) theta2(t) r1 r2 m1 m2 g 

% Initalize Kinematic Constraints
x1 = r1*sin(theta1);
y1 = -r1*cos(theta1);
x2 = x1+r2*sin(theta2+theta1);
y2 = y1-r2*cos(theta2+theta1);

% Initalize Velocity Equations
x1dot=diff(x1);
y1dot=diff(y1);
x2dot=diff(x2);
y2dot=diff(y2);

% Solve for Potential Energy
pe=m1*g*y1+m2*g*y2;

%Solve for Kinetic Energy
ke=1/2*m1*(x1dot^2+y1dot^2)+1/2*m2*(x2dot^2+y2dot^2);
simplify(ke);

% Solve for L=KE-PE
l=ke-pe;
simplify(l);

% Initalize Variables
syms theta1dot theta2dot

% Solve for partial derivitive of L by theta1dot
eqn1=subs(diff(subs(l(t), diff(theta1(t)), theta1dot), theta1dot), theta1dot, diff(theta1(t)));
eqn1=simplify(eqn1);

% Solve for the time derivative of eqn1 
eqn1=diff(eqn1);
eqn1=simplify(eqn1);

% Solve for partial derivitive of L by theta1
a=simplify(diff(l,theta1));

% Solve for first equation of motion
eqn1=simplify(eqn1-a)

% Solve for partial derivitive of L by theta2dot
eqn2=subs(diff(subs(l(t), diff(theta2(t)), theta2dot), theta2dot), theta2dot, diff(theta2(t)));
eqn2=simplify(eqn2);

% Solve for the time derivative of eqn2 
eqn2=diff(eqn2);
eqn2=simplify(eqn2);

% Solve for partial derivitive of L by theta2
b=simplify(diff(l,theta2));

% Solve for second equation of motion
eqn2=simplify(eqn2-b)

r1 = 1;
r2 = 1.5;
m1 = 2;
m2 = 1;
g = 9.8;

[V,S] = odeToVectorField(eqn1,eqn2)

M = matlabFunction(V,'vars',{'t','Y'})

initCond = [pi/4 0 pi/6 0];
sols = ode45(M,[0 10],initCond);

plot(sols.x,sols.y)
legend('\theta_2','d\theta_2/dt','\theta_1','d\theta_1/dt')
title('Solutions of State Variables')
xlabel('Time (s)')
ylabel('Solutions (rad or rad/s)')

Answer 1

[Walter Roberson]

If you define a symbolic variable and use it in an expression, and later assign a numeric value to the same name, then that does not change the symbolic variable that is in the expression. The situation is much like how if you

A = 3

B = A*2 + 1

A = 9

then B does not suddenly change from 3*2+1 to instead evaluate to 9*2+1 . The value of A as of the time that B was defined is what is used in B, and changes to A after that do not affect B. The same is true for if it had been syms A instead of assigning a numeric value to A -- changing A afterwards would not change the place the symbolic version was already used.

% Initalize Variables
syms theta1(t) theta2(t) r1 r2 m1 m2 g 

% Initalize Kinematic Constraints
x1 = r1*sin(theta1);
y1 = -r1*cos(theta1);
x2 = x1+r2*sin(theta2+theta1);
y2 = y1-r2*cos(theta2+theta1);

% Initalize Velocity Equations
x1dot=diff(x1);
y1dot=diff(y1);
x2dot=diff(x2);
y2dot=diff(y2);

% Solve for Potential Energy
pe=m1*g*y1+m2*g*y2;

%Solve for Kinetic Energy
ke=1/2*m1*(x1dot^2+y1dot^2)+1/2*m2*(x2dot^2+y2dot^2);
simplify(ke);

% Solve for L=KE-PE
l=ke-pe;
simplify(l);

% Initalize Variables
syms theta1dot theta2dot

% Solve for partial derivitive of L by theta1dot
eqn1=subs(diff(subs(l(t), diff(theta1(t)), theta1dot), theta1dot), theta1dot, diff(theta1(t)));
eqn1=simplify(eqn1);

% Solve for the time derivative of eqn1 
eqn1=diff(eqn1);
eqn1=simplify(eqn1);

% Solve for partial derivitive of L by theta1
a=simplify(diff(l,theta1));

% Solve for first equation of motion
eqn1=simplify(eqn1-a)

eqn1(t) = 

% Solve for partial derivitive of L by theta2dot
eqn2=subs(diff(subs(l(t), diff(theta2(t)), theta2dot), theta2dot), theta2dot, diff(theta2(t)));
eqn2=simplify(eqn2);

% Solve for the time derivative of eqn2 
eqn2=diff(eqn2);
eqn2=simplify(eqn2);

% Solve for partial derivitive of L by theta2
b=simplify(diff(l,theta2));

% Solve for second equation of motion
eqn2=simplify(eqn2-b)

eqn2(t) = 

r1 = 1;
r2 = 1.5;
m1 = 2;
m2 = 1;
g = 9.8;

[V,S] = odeToVectorField(subs(eqn1),subs(eqn2))

V = 

S = 

M = matlabFunction(V,'vars',{'t','Y'})

M = function_handle with value:
    @(t,Y)[Y(2);(sin(Y(3)).*-2.94e+2+sin(Y(1)+Y(3)).*1.96e+2+sin(Y(1)).*Y(2).^2.*1.5e+1+sin(Y(1)).*Y(4).^2.*3.5e+1-cos(Y(1)).*sin(Y(3)).*1.96e+2+cos(Y(1)).*sin(Y(1)+Y(3)).*9.8e+1+cos(Y(1)).*sin(Y(1)).*Y(2).^2.*1.0e+1+cos(Y(1)).*sin(Y(1)).*Y(4).^2.*2.0e+1+sin(Y(1)).*Y(2).*Y(4).*3.0e+1+cos(Y(1)).*sin(Y(1)).*Y(2).*Y(4).*2.0e+1)./(cos(Y(1)).^2.*1.0e+1-3.0e+1);Y(4);((sin(Y(3)).*-2.94e+2+sin(Y(1)).*Y(2).^2.*1.5e+1+sin(Y(1)).*Y(4).^2.*1.5e+1+cos(Y(1)).*sin(Y(1)+Y(3)).*9.8e+1+cos(Y(1)).*sin(Y(1)).*Y(4).^2.*1.0e+1+sin(Y(1)).*Y(2).*Y(4).*3.0e+1).*(-1.0./1.0e+1))./(cos(Y(1)).^2-3.0)]

initCond = [pi/4 0 pi/6 0];
sols = ode45(M,[0 10],initCond);

plot(sols.x,sols.y)
legend('\theta_2','d\theta_2/dt','\theta_1','d\theta_1/dt')
title('Solutions of State Variables')
xlabel('Time (s)')
ylabel('Solutions (rad or rad/s)')