We propose a numerical method (implemented in Matlab) for computing approximate Fekete points on compact multivariate domains. It relies on the search of maximum volume submatrices of Vandermonde matrices computed on suitable discretization meshes, and uses a simple greedy algorithm based on QR factorization with column pivoting. The method gives also automatically an algebraic cubature formula, provided that the moments of the underlying polynomial basis are known. Numerical tests are presented for the interval and the square, which show that approximate Fekete points are well suited for polynomial interpolation and cubature.