Unfortunately, there isn't a painless way to do this. That is, it's likely that writing new code to get eigenvectors given the eigenvalues will take you longer to get right, then just re-running eigs.
If you had the eigenvectors and needed the eigenvalues, that's quite painlessly possible using x'*A(x) for each eigenvector x.
The other way around, the usual way to do this would be to solve (A - lambda_mod*I) * x = rand, where lambda_mod is slightly different from lambda to avoid a singular linear system. You might need to do this in multiple iterations (inverse iteration as Torsten mentioned). But with your matrix being given as a function handle, I imagine solving a linear system with it is non-trivial.
