Isolated pitching and plunging airfoils have been extensively studied in literature. The effects of proximity of sections of flapping wings to each other and to solid walls have received much less attention. The optimal spacing and program of motion of flapping foils can generate thrust instead of drug. Thrust generation studies of flapping wings in proximity to other wing(s) and/or walls are needed for micro air vehicles, marine propulsion and generation of hydroelectric power. A numerical approach must handle complex multi--body dynamics while avoiding use of body-fitted grids needed for traditional finite-volume computations, which would require domain re-meshing at all time steps. The adopted and tested in the current study computational approach is based on particle-based kinetic Lattice Boltzmann method (LBM). Commercial LBM-based CFD software XFlow has been adopted, validated and tested. The local level of refinement of the underlying lattice is selected with one parameter controlling the refinement in near-surface boundary layer and another parameter controlling refinement in street of vortices behind the airfoil. Governing parameters of single and paired flapping foils investigated include profile of foil(s) and thickness, amplitude of plunging, angular amplitude of pitching, phase difference between plunging and pitching, proximity to rigid wall and distance between foils in tandem.