Wave propagation in elastic metamaterials has an extremely low energy dissipation [1], which is a key benefit for passive (or near to zero) neuromorphic information processors. Such complicated metamaterials can be designed using neural networks [2,3]. However, these neural networks are too complex because they have a huge number of neurons which contain nonlinear resonators. Therefore, it is not reasonable to use FEM. Hence, we need to develop a new efficient numerical method to model the nonlinear response of elastic metamaterials. In this research, we developed a new model order reduction technique which is based on Wannier functions (inspired from Density Functional Theory). The two main advantages of this method are the sparsity of the reduced nonlinear tensors and the localization of the base vectors. The latter advantage makes the oscillation of each oscillator depending on itself and its neighborhood. Hence, instead of considering the entire lattice, simulating the small clusters of sites is sufficient. By using the basis of maximally-localized functions we can simulate highly complex metamaterials efficiently. Additionally, by adding a quadratic manifold technique [4,5] higher-order modes are statically condensed, even without increasing the computational costs significantly. The method is applied to both dynamics and static problems. Non-trivial properties such as buckling can be captured. We implemented the method on designing different metamaterials and comparing the results with high-fidelity simulations.