The neural-based approaches inspired by biological neural mechanisms of locomotion are becoming increasingly popular in robot control. This paper investigates a systematic method to formulate a Central Pattern Generator (CPG) based control model for multimodal swimming of a multi-articulated robotic fish with flexible pectoral fins. A CPG network is created to yield diverse swimming in three dimensions by coupling a set of nonlinear neural oscillators using nearest-neighbor interactions. In particular, a sensitivity analysis of characteristic parameters and a stability proof of the CPG network are given. Through the coordinated control of the joint CPG, caudal fin CPG, and pectoral fin CPG, a diversity of swimming modes are defined and successfully implemented. The latest results obtained demonstrate the effectiveness of the proposed method. It is also confirmed that the CPG-based swimming control exhibits better dynamic invariability in preserving rhythm than the conventional body wave method. Robotic fish, inspired by fish in nature, have received increasing attention in recent years because of potential aquatic-related applications such as exploration, surveillance, transportation, and mobile sensing [1][2][3][4][5][6][7]. Many scientists and engineers have carried out studies on robotic fish, focusing not only on the hydrodynamics mechanism and bio-inspired practice of fishlike swimming, but also on motion planning, control, and optimization of the fishlike robots [2]. The current swimming control methods, from the perspective of cybernetics, tend to fall into two primary categories: bio-inspired and non-bio-inspired (conventional). The former is nourished by an abundance of biological knowledge of fish or other animals, while the latter relates to deriving control laws from a combined analysis of multi-body dynamics and kinematics. To achieve flexible swimming control, conventional locomotion control usually seeks to calculate propulsive forces and to determine the accompanying movements. Once the dynamic equation attached to each propulsive component has been identified via synthesizing all involved dynamic equations, the control laws can be obtained in terms of the resultant forces and kinematic parameters of the moving joints. As a rule, oscillatory frequency, amplitude, as well as phase difference between adjacent joints are usually extracted as the locomotion-related characteristics significantly affecting swimming performance. Many robotic fishes are controlled in such a manner, for example, the well-known Robotuna [8], the pectoral-fin-driven robotic Blackbass [9], and the lifelike Mitsubishi robotic fish [10]. However, two issues arise in determining appropriate swimming parameters: (1) the difficulty in accurately modeling the hydrodynamic interaction between the oscillating fish body and its surrounding fluid, and (2) an infinity of possible solutions in the hyper-redundant planar kinematic