if(mod(t,1000)==0)
figure(1), clf, box
    imagesc([min(min(X)) max(max(X))],[min(min(Y)) max(max(Y))], T(iny,inx)), axis equal, ...
        axis tight, drawnow, colorbar
        ylabel('y [m]')
    xlabel('x [m]'), pause(0.1)

figure(2), clf, box
    dTdywall = (T(1,inx)-T(2,inx))/dy; 
    Tmean = sum(U.*T(iny,inx))./sum(U);
    qwall = k*dTdywall; 
    plot(X(2,:), qwall,'k-','Linewidth', 2)
    xlabel('x')
    ylabel('Wall heat flux [W/m^2]')

figure(3), clf, box
    plot(X(2,:), Ly*2*dTdywall./(Ttop-Tmean),'k-','Linewidth', 2)
    xlabel('x')
    ylabel('Nusselt number')
    
figure(4), clf, box
    plot(X(2,:), Tmean,'k-','Linewidth', 2)
    xlabel('x')
    ylabel('Temperature [K]')
     
spreadfigures

end