We study the response of topological insulator films to strong magnetic and electric fields beyond the linear response theory. As a model, we use a three-dimensional lattice Wilson-Dirac Hamiltonian where we simultaneously introduce both magnetic field through the Peierls substitution and electric field as a potential energy depending on lattice coordinate. We compute the electron energy spectrum by numerically diagonalizing this Hamiltonian and obtain quantized magnetoelectric polarizability. In addition, we find that the magnetoelectric effect vanishes as the film width decreases due to the hybridization of surface wave functions. Furthermore, applying a gate voltage between the surfaces, we observe several quantized plateaus of θ term, which are mainly determined by the Landau level structures on the top and bottom surfaces.