The problem here is that you have defined 4 equations and trying to solve for 5 unknowns (beta1, beta2 beta3 M2, M3) that is incompatible. Thus, one of the unknowns to be known/defined. Here is corrected code with this alteration.
syms M3 beta1 beta2 beta3
M1 = 2.4;
M2 = 4.2; % Defined and M3 remains as unknown
M4 = 1.3;
Gamma = 1.4;
Gamma_P = Gamma+1;
Gamma_M = Gamma-1;
%Mach #'s upstream of normal shock
Mn1 = M1*sind(beta1);
Mn2 = M2*sind(beta2);
Mn3 = M3*sind(beta3);
%Deflection Angles
theta1 = atand(2*(cotd(beta1))*(((M1^2)*(sind(beta1))^2)-1)/((2+((M1^2)*(Gamma+cosd(2*beta1))))));
theta2 = atand(2*(cotd(beta2))*(((M2^2)*(sind(beta2))^2)-1)/((2+((M2^2)*(Gamma+cosd(2*beta2))))));
theta3 = atand(2*(cotd(beta3))*(((M1^2)*(sind(beta3))^2)-1)/((2+((M3^2)*(Gamma+cosd(2*beta3))))));
% Mach #'s
M2x = ((sqrt(((Mn1^2*Gamma_M)+2)/((2*Gamma*Mn1^2)-Gamma_M))))/(sind((beta1)-(theta1)));
M3x = ((sqrt(((Mn2^2*Gamma_M)+2)/((2*Gamma*Mn2^2)-Gamma_M))))/(sind((beta2)-(theta2)));
M4x = ((sqrt(((Mn3^2*Gamma_M)+2)/((2*Gamma*Mn3^2)-Gamma_M))))/(sind((beta3)-(theta3)));
solve = vpasolve([Mn1==Mn2, Mn1==Mn3, M2==M2x, M3==M3x],[beta1 beta2 beta3 M3],[10 20 17 3.5])
