syms Ta Tc alpha gammar Tz xiz QH deltaQH M1 betaz
mcv=19.806+(0.00419/2).*Ta
mcphayphayvz=19.806+(0.00419/2).*Tz+0.5.*Tz.*((360.34+252.4*alpha)^(10^-5))
mcphayvc=(mcv+gammar.*mcphayphayvz)./(1+gammar)
eqn = Tz==(((xiz*(QH-deltaQH))./(M1.*(1+gammar)))+mcphayvc.*Tc)./(betaz.*mcphayphayvz);
Tzsol = solve(eqn, Tz)
Note that Ta Tc gammar must be scalar values for this purpose. If they were arrays then the solve() would be trying to find one single Tz that solved all of the values simultaneously.
If you are trying to get a solution for each corresponding Ta, Tc, gammar, then you cansubs(Tzsol, {Ta, Tc, gammar}, {Ta_array, Tc_array, gammar_array})
Also note that a pair of solutions are found. It looks like it might be two solutions for a quadratic.