Studies on hydrodynamic characteristics of viscous incompressible flows around flexible hydrofoils are of practical importance for the design and performance optimization of marine structures such as ship rudders and stabilizing fins. The aim of this paper is to extend a radial basis function based ghost cell method to simulate flows around single or multiple flexible moving hydrofoils in array arrangement. The numerical model is based on a ghost cell finite difference method for considering the influence of the immersed boundaries on the flow. Also, a compact supported radial basis function (CSRBF) is introduced to track the complex flexible boundary with some controlling points of the body surface. Based on the present method, the uniform flow around a flexible hydrofoil swimming like a fish is simulated. Good grid convergence of drag and lift coefficients demonstrates the accuracy and reliability of the present method. Also, the hydrodynamics patterns of the flexible hydrofoil under different oscillation frequencies are studied. Further, the thrust generation mechanism of the hydrofoil is explained. Afterwards, flows around the undulating hydrofoils in array arrangement are simulated. The force coefficients and wake patterns under different distances and oscillation frequencies are investigated. It is observed that the thrust coefficients of the hydrofoils under narrow arrangement and high oscillation frequencies have significant amplification effects. In addition, the critical frequency at zero thrust reduces.