disp(strcat('Maximum Courant number (Co=dt*U/dx): ', num2str(dt*max(max(U))/dxmin)))
disp(strcat('Maximum CFL number (CFL=dt*alpha/dx^2): ', num2str(dt*alpha/dxmin^2)))

Re=mean(U(:,1))*Ly/alpha;