If I understand correctly, the column used to get the x-coordinates out of the matrix lnRe_phi, needs to be 1, 2, 3, 4, 5 when i5 is 0, 4, 8, 12, 16. If that's right, then something like this would work:
% 26-by-5 matrix of x coordinates
lnRe_phi = linspace(0,5,26).'+(0:4)
% 26-by-20 matrix of y coordinates
lnNu_r = randn(26,20);
N1 = 4;
for i5 = 0:N1:(N1*N1)
%h1 = h(:,i5+1);
%h2 = h(:,i5+2);
%h3 = h(:,i5+3);
%h4 = h(:,i5+4);
lnNu_r1 = lnNu_r(:,i5+1);
lnNu_r2 = lnNu_r(:,i5+2);
lnNu_r3 = lnNu_r(:,i5+3);
lnNu_r4 = lnNu_r(:,i5+4);
x_column = i5/N1+1;
% lnNu vs lnRe for effect of Ctflow
figure(i5+2)
plot(lnRe_phi(:,x_column),lnNu_r1,'+')
hold on
plot(lnRe_phi(:,x_column),lnNu_r2,'+')
plot(lnRe_phi(:,x_column),lnNu_r3,'+')
plot(lnRe_phi(:,x_column),lnNu_r4,'+')
hold off
title('Variation of ln(Local Nusselt Number) vs ln(Local Reynolds Number)')
xlabel('ln(Re_\phi)')
ylabel('ln(Nu_r)')
legend('{C}_{{w}_{1}}','{C}_{{w}_{2}}','{C}_{{w}_{3}}','{C}_{{w}_{4}}')
grid
end
