how to solve BdG Equations numerically

[Jason Wang]

For example,
how can one get the BdG equations for Bose Hubbard model
Then how to solve it numerically?

Anyone may help?