[en] A method is described to compute invariant tori in phase space for calssical non-integrable Hamiltonian systems. Our procedure is to solve the Hamilton-Jacobi equation stated as a system of equations for Fourier coefficients of the generating function. The system is truncated to a finite number of Fourier modes and solved numerically by Newton's method. The resulting canonical transformation serves to reduce greatly the non-integrable part of the Hamiltonian. In examples studied to date the convergence properties of the method are excellent, even near chaotic regions and on the separatrices of isolated broad resonances. We propose a criterion for breakup of invariant surfaces, namely the vanishing of the Jacobian of the canonical transformation to new angle variables. By comparison with results from tracking, we find in an example with two nearly overlapping resonances that this criterion can be implemented with sufficient accuracy to determine critical parameters for the breakup ('transition to chaos') to an accuracy of 5 to 10%