It looks like you are trying to simplify the equation by substituting trigonometric identities for the cosines and sines of the angles theta1, theta2, and theta3. You need to substitute these trigonometric identities into the expression before expanding and simplifying. To do this, you can use the subs function again, as in the following example:
eqT_01 = subs(eqT_01, [c1 s1 c2 s2 c3 s3 c23 s23], [cos(theta1) sin(theta1) cos(theta2) sin(theta2) cos(theta3) sin(theta3) cos(theta2)*cos(theta3)-sin(theta2)*sin(theta3) cos(theta2)*sin(theta3)+sin(theta2)*cos(theta3)])
Once you have performed this substitution, you can then expand and simplify the equation using the expand and simplify functions, as you have done already.
Here is the complete code with the changes I suggested:
syms s [1 6]; syms c [1 6]; syms a2 a3 d3 d4; syms R_1 [1 3];syms R_2 [1 3];syms R_3 [1 3] ;
syms Px Py Pz;
syms theta1 theta2 theta3 theta4 theta5 theta6;
syms c23 s23;
T_01_14 = a2*c2 - Py*s1 - Px*c1 + a3*c23 - d4*s23
T_01_34 = -Pz - c23*d4 - a2*s2 - a3*s23
T_01_24 = Py*c1 - d3 - Px*s1
eqT_01 = T_01_14^2 + T_01_34^2 + T_01_24^2
% Substitute trigonometric identities for the cosines and sines of the angles
eqT_01 = subs(eqT_01, [c1 s1 c2 s2 c3 s3 c23 s23], [cos(theta1) sin(theta1) cos(theta2) sin(theta2) cos(theta3) sin(theta3) cos(theta2)*cos(theta3)-sin(theta2)*sin(theta3) cos(theta2)*sin(theta3)+sin(theta2)*cos(theta3)])
% Expand and simplify the equation
eqT_01 = expand(eqT_01)
eqT_01 = simplify(eqT_01)
