The stability of the immersed boundary (IB) method is a challenge in simulating fluid-structure interaction (FSI) problems, where time step constraints are significantly stricter than in pure fluid simulations. We propose a novel unconditionally stable scheme for the immersed boundary finite element (IBFE) method. The structure is handled implicitly and characterized by strain energy functions, rather than being modeled as fibers or membranes. Through energy estimate, we prove the unconditional stability of the fully discrete approximation in the absence of the convective term. In real simulations of cardiac mechanics problems, the time step is much larger, only limited by the Courant–Friedrichs–Lewy (CFL) condition of the fluid. The novelty of this work lies in the combination of dual interpolation and distribution operators, the Jacobian-free Newton-Krylov (JFNK) method for solving nonlinear algebraic systems, and the semi-Lagrangian method for handling the convective term. To validate the effectiveness and accuracy of our approach, we present various benchmarks and conduct a quasi-static simulation of a three-dimensional real left ventricular model. We have shown that the numerical stability of our scheme is very robust even with much larger time step compared to conventional explicit IB methods. Our work paves the way for further works on efficient solvers of the IBFE method.