I'm sorry I don't have time to fully understand your code, from a quick look I think your problem may be in your function at line 98
y2_2nd = @(t,q,dqdt) (Q_q-C_q.*dqdt.^2)./J_q;
along with the later line 102
k_12 = y2_2nd(t(i),q(i),dqdt(i));
Your function contains the terms Q_q and C_q both of which are vectors. This means that y2_2nd returns a vector and as result k_12 is also a vector. I'm thinking that you want k_12 to be a scalar constant.
Similar problems then occur in the later terms lines 103-108 so that finally at line 110 you are trying to assign a scalar, q(i+1) to a vector, since, k_21, k_31, and k_41 are all vectors so the right hand side of line 110 evaluates to a vector.