You create T_g T_b T_c as numeric variables.
That is your only symbolic variable, and it is a scalar.
Eq1=(alpha_g*I(i))+(h_ge*(T_e(i)-T_g(i)))+(U_gc*(T_c(i)-T_g(i)));
Eq2=Fc*alpha_c*tha_g*I(i)*((1/Fc)-etha_ref-(etha_ref*beta*T_ref)+(etha_ref*beta*T_c(i)))+(U_gc*(T_g(i)-T_c(i)))+(U_cb*(T_b(i)-T_c(i)));
Eq3=(((Fc*tha_c)-Fc+1)*alpha_b*tha_g*I(i))+(U_cb*(T_c(i)-T_b(i)))+(U_bp*(T_p-T_b(i)));
The first time through the loop (at least), T_g(i), T_b(i), T_c(i) are all numeric zeros.
Eq1 is independent of the symbolic variable T_p so it is effectively a constant the first time through. The same is true about Eq2, independent of all symbolic variables (at least in the first pass). Eq3 does depend upon T_p, in a relatively trivial way.
[T_g,T_c,T_b]=solve([Eq1,Eq2,Eq3],[T_g,T_c,T_b]);
At the time of the first pass through the loop, T_g, T_c, T_b are all numeric vectors of zeros, but they appear in the slot reserved for lists of variables to solve for. You have three equations in one unknown, and two of the equations are numeric constants. If the constants just happen to equal 0 exactly then it might be possible to solve for T_p... but not while the [T_g,T_c,T_b] are numeric.
In MATLAB, if you are not using inequalities you should have the same number of equations as you have variables being solved for.
If, hypothetically, your solve() worked, then you overwrite all of T_g, T_c, and T_b . Unless they happen to get overwritten with at least N values, you would have trouble in the next iteration of the loop when you try to use T_g(i) and so on.