The axisymmetric squeeze flow of viscoplastic materials is examined using either the original Bingham constitutive equation or the approximate model suggested by Papanastasiou. Previous theoretical analyses of this problem, using the standard lubrication approximation, have led to conflicting results, whereby the material around the plane of symmetry must both behave as unyielded solid and translate radially with a nonuniform velocity. In these analyses, the ratio of half the distance between the plates to their radius, epsilon, is taken to be small. With a qualitative analysis we show that unyielded material must exist only around the two stagnation points of flow at the center of the disks and must cover only a fraction of the axis of symmetry. It is also shown that normal forces must be included around the plane of symmetry in order to resolve the earlier paradox of approximate analyses. Our converged numerical results verify these ideas and demonstrate further that a typical lubrication analysis is bound to fail in this problem, since the velocity and pressure fields remain two-dimensional, even when epsilon much less than 1. Using the original Bingham model we compute simultaneously the shape of the yield surface, the velocity and pressure fields employing the Galerkin/finite element methodology. Using the Papanastasiou model computations become simpler avoiding the a priori calculation of the yield surface. However, if its exponential parameter is taken to be sufficiently large, the computed velocity and pressure fields are in very good agreement with those obtained from the original Bingham model. The results are primarily affected by the Bingham number that measures the magnitude of the yield stress with respect to the viscous stresses. As this number increases, large departures from the corresponding Newtonian solution are obtained and limited flow and deformation of the material is predicted.