The classical n-vector ϕ{4} model with O(n) symmetrical Hamiltonian H is considered in a ∞{2}×L slab geometry bounded by a pair of parallel free surface planes at separation L. Standard quadratic boundary terms implying Robin boundary conditions are included in H. The temperature-dependent scaling functions of the excess free energy and the thermodynamic Casimir force are computed in the large-n limit for temperatures T at, above, and below the bulk critical temperature T_{c}. Their n=∞ limits can be expressed exactly in terms of the spectrum and eigenfunctions of a self-consistent one-dimensional Schrödinger equation. This equation is solved by numerical means for two distinct discretized versions of the model: in the first ("model A"), only the coordinate z across the slab is discretized and the integrations over momenta conjugate to the lateral coordinates are regularized dimensionally; in the second ("model B"), a simple cubic lattice with periodic boundary conditions along the lateral directions is used. Renormalization-group ideas are invoked to show that, in addition to corrections to scaling ∝L{-1}, anomalous ones ∝L{-1}lnL should occur. They can be considerably decreased by taking an appropriate g→∞ (T_{c}→∞) limit of the ϕ{4} interaction constant g. Depending on the model A or B, they can be absorbed completely or to a large extent in an effective thickness L_{eff}=L+δL. Excellent data collapses and consistent high-precision results for both models are obtained. The approach to the low-temperature Goldstone values of the scaling functions is shown to involve logarithmic anomalies. The scaling functions exhibit all qualitative features seen in experiments on the thinning of wetting layers of {4}He and Monte Carlo simulations of XY models, including a pronounced minimum of the Casimir force below T_{c}. The results are in conformity with various analytically known exact properties of the scaling functions.