Computation of Maximally Localized Wannier Functions using a simultaneous diagonalization algorithm