Flexible intelligent metasurfaces (FIMs) show great potential for improving the wireless network capacity in an energy-efficient manner. An FIM is a soft array consisting of several low-cost radiating elements. Each element can independently emit electromagnetic signals, while flexibly adjusting its position even perpendicularly to the overall surface to `morph' its 3D shape. More explicitly, compared to a conventional rigid antenna array, an FIM is capable of finding an optimal 3D surface shape that provides improved signal quality. In this paper, we study point-to-point multiple-input multiple-output (MIMO) communications between a pair of FIMs. In order to characterize the capacity limits of FIM-aided MIMO transmissions over frequency-flat fading channels, we formulate a transmit optimization problem for maximizing the MIMO channel capacity by jointly optimizing the 3D surface shapes of the transmitting and receiving FIMs as well as the MIMO transmit covariance matrix, subject to the total transmit power constraint and to the maximum perpendicular morphing range of the FIM. To solve this problem, we develop an efficient block coordinate descent (BCD) algorithm. The BCD algorithm iteratively updates the 3D surface shapes of the FIMs and the transmit covariance matrix, while keeping the other fixed, to find a locally optimal solution. Numerical results verify that FIMs can achieve higher MIMO capacity than that of the conventional rigid arrays. In particular, the MIMO channel capacity can be doubled by the proposed BCD algorithm under some setups.