Since thetac is an array of 15 elements, you are trying to solve a system of 15 incompatible equations, therefore, there is no answer. You should solve the equation for each value of thetac.
syms theta1 theta2 gama M l
eqn = theta1 == 2*atan((2*tan(theta2/2)*tan(gama/2 - theta2/2) + (4*tan(theta2/2)^2 - 4*M^2*tan(theta2/2)^2 - 4*M^2*tan(gama/2 - theta2/2)^2 - M^2*gama^2*l^2 - 8*M^2*tan(theta2/2)*tan(gama/2 - theta2/2) - 4*M^2*gama*l*tan(theta2/2) - 4*M^2*gama*l*tan(gama/2 - theta2/2) - M^2*gama^2*l^2*tan(theta2/2)^2*tan(gama/2 - theta2/2)^2 + 2*M^2*gama^2*l^2*tan(theta2/2)*tan(gama/2 - theta2/2) + 4*M^2*gama*l*tan(theta2/2)*tan(gama/2 - theta2/2)^2 + 4*M^2*gama*l*tan(theta2/2)^2*tan(gama/2 - theta2/2) + 4)^(1/2) - 2)/(2*tan(theta2/2) + 2*M*tan(gama/2 - theta2/2) + 2*M*tan(theta2/2) + M*gama*l - M*gama*l*tan(theta2/2)*tan(gama/2 - theta2/2)) - (2*tan(theta2/2)*tan(gama/2 - theta2/2))/(2*tan(theta2/2) + 2*M*tan(gama/2 - theta2/2) + 2*M*tan(theta2/2) + M*gama*l - M*gama*l*tan(theta2/2)*tan(gama/2 - theta2/2)));
F=1.06;
gama=pi/F;
thetac=0:0.2:gama;
l=0.2181818;
M=0.960615;
for k=1:numel(thetac)
theta2=thetac(k);
seqn = subs(eqn);
sol(k) = vpasolve(seqn)
end
alpha_1=sol./thetac;
plot(thetac,alpha_1)
