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 Tc. 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−1lnL should occur. They can be considerably decreased by taking an appropriate g→∞ (Tc→∞) 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 Leff=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 4He and Monte Carlo simulations of XY models, including a pronounced minimum of the Casimir force below Tc. The results are in conformity with various analytically known exact properties of the scaling functions.