Abstract. The Immersed Boundary Method (IB) is known as a powerful technique for the numerical solution of fluid-structure interaction problems as, for instance, the motion and deformation of viscoelastic bodies immersed in an external flow. It is based on the treatment of the flow equations within an Eulerian framework and of the equations of motion of the immersed bodies with respect to a Lagrangian coordinate system including interaction equations providing the transfer between both frames. The classical IB uses finite differences, but the IBM can be set up within a finite element approach in the spatial variables as well (FE-IB). The discretization in time usually relies on the Backward Euler (BE) method for the semidiscretized flow equations and the Forward Euler (FE) method for the equations of motion of the immersed bodies. The BE/FE FE-IB is subject to a CFL-type condition, whereas the fully implicit BE/BE FE-IB is unconditionally stable. The latter one can be solved numerically by Newton-type methods whose convergence properties are dictated by an appropriate choice of the time step size, in particular, if one is faced with sudden changes in the total energy of the system. In this paper, taking advantage of the well developed affine covariant convergence theory for Newton-type methods, we study a predictor-corrector continuation strategy in time with an adaptive choice of the continuation steplength. The feasibility of the approach and its superiority to BE/FE FE-IB is illustrated by a representative numerical example.Key words. finite element immersed boundary method, fully implicit scheme, predictorcorrector continuation, red blood cells AMS subject classifications. 65H20, 65M60, 74L15, 76D05, 92C101. Introduction. A computationally attractive methodology for the numerical simulation of the motion and deformation of elastic and viscoelastic bodies in external flows is the Immersed Boundary Method (IB), which has been originally developed by Peskin [28] and further studied in [11,29,30,31,33]. The IBM uses an Eulerian coordinate system for the flow equations and Lagrangian coordinates for the boundary of the immersed bodies together with appropriate interaction equations to transform Eulerian to Lagrangian quantities and vice versa. The interaction equations feature multidimensional Dirac delta functions that have to be approximated appropriately within a finite difference approach. More recently, a variational formulation of the IBM has been provided in [8] and [9, 10] as a basis for a finite element realization referred to as the Finite Element Immersed Boundary Method (FE-IB). Both for the classical IB and the FE-IB, the most common approach with regard to discretization in time is to use the Backward Euler (BE) method for the flow equations and the Forward Euler (FE) method for the equation describing the motion and deformation of the immersed bodies which gives rise to the BE/FE IB and BE/FE FE-IB, respectively. However, these schemes typically require a CFL-type condition (cf., e.g., [9,16]). B...