1 - This is most certainly due to the choice of boundary condition. The code uses bounce-back for no-slip conditions and Zou/He for the lid, which, I think, is not the numerically most stable combination. It was chosen for the sake of illustration, as both boundary conditions are commonly used. You can find other conditions in the literature, though, as for example here.
2 - You are completely right. Choosing a value of 400’000 for the total number of iterations is arbitrary, and there is no guarantee that the simulation has converged to a stationary state. Actually, the purpose of this Matlab script is to illustrate how to implement a simple LB code, and not to provide a full study of a physical problem. If you really want to do physics with this code (but personally, I wouldn’t, as it is quite slow; why not use OpenLB instead?) you’ll probably want to replace the for loop by a loop which breaks when a stationary state is reached.
3 - The simulation runs in lattice units. You can easily convert to any physical units, as explained here.