Not entirely clear to me what you want, but perhaps something like the following, where the edge equations are taken outside the loop:
% Arbitrary data
Qc = 0;
Qh = 0;
q0 = 1000;
Pc = 1;
Ph = 1;
qw = 1;
Rc = 1;
R1 = 1;
Rg0 = 1;
Rh0 = 1;
Nc = 1;
T0 = 10;
Ntot = 10;
%unknowns
Qm = zeros(Ntot,1);
T = zeros(Ntot,1);
%system 1
Qm(1) = -Qc + q0 - Pc - qw/2;
T(1) = T0 - Rc*Qc;
for ii=2:Ntot-1
%system i
T(ii) = T(ii-1) + R1/Nc*q0 - R1/Nc*Qm(ii-1);
Qm(ii) = Qm(ii-1) + (T0 - T(ii))*Nc/Rg0;
end
%system N
Qm(Ntot) = Qh - (Ph - q0 + qw/2);
T(Ntot) = T0 + Rh0*Qh;
plot(Qm,T),grid
