Consider first the situation in 1D: You have N points for with conditions from the eigenvalue equation and two boundary conditions, which makes N+2 conditions in total. Yet, there are only N degrees of freedom, which are given by the expansion coefficents of the collocation functions. This makes an overdetermined system which can't be solved in general.
In order to make the system solvable, there are basically two strategies: either drop two conditions completely (--usually not the BCs), or use some kind of weighting scheme such as least squares. The former satisfies the N selected constraints exactly, while the latter tries to satisfy all constraints in an approximate fashion.
You asked for the former scheme. For this, one faces the problem how to distribute the N constraints to the N points. Dirichlet condition at a given collocation point are handled quite naturally by simply fixing the value at this point, even more so for zero Dirichlet conditions (which allow to completely remove the corresponding collocation function from the basis expansion).
For more general conditions, e.g. Dirichlet between points, Neumann or Robin conditions, The selection of the degree of freedom which should handle the BC is not that straightforward. Thus, you can either choose aby point (usually the boundary point), or you apply a more clever scheme such as rectangular spectral collocation by Driscoll and Hale. The latter basically applies a projection down to the space of polynomial function of degree N-2, wgich gives the freedom to add the two BCs without runnimg into an overdetermined system. Rectangular spectral collocation also was shown to behave more consistently and better conditioned. I'd suggest you have a look at the original paper and try to apply it to your problem. Nothing more concrete, sorry, as it will be a bit of fiddling.
EDIT: the problem will then boil down to mapping the cartesian product grid of total size N down to a grid of N-2. That's conceptually easy in 1D, but in 3D it's again a cpmplex step...