Research on fish locomotion has made extensive progress towards a better understanding of how fish control their flexible body and fin for propulsion and maneuvering. Although the biologically flexible fish fins are believed to be one of the most important features to achieve optimal swimming performance, due to the limitations of the existing numerical modeling tool, studies on a deformable fin with a non-uniformly distributed stiffness are rare. In this work, we present a fully coupled fluid–structure interaction solver which can cope with the dynamic interplay between flexible aquatic animal and the ambient medium. In this tool, the fluid is resolved by solving Navier–Stokes equations based on the finite volume method with a multi-block grid system. The solid dynamics is solved by a nonlinear finite element method. A sophisticated improved IQN-ILS coupling algorithm is employed to stabilize solution and accelerate convergence. To demonstrate the capability of the developed Fluid–Structure-Interaction solver, we investigated the effect of five different stiffness distributions on the propulsive performance of a caudal peduncle-fin model. It is shown that with a non-uniformly distributed stiffness along the surface of the caudal fin, we are able to replicate similar real fish fin deformation. Consistent with the experimental observations, our numerical results also indicate that the fin with a cupping stiffness profile generates the largest thrust and efficiency whereas a heterocercal flexible fin yields the least propulsion performance but has the best maneuverability.