Integral formulations can be convenient for computing eddy currents in complicated electromagnetic systems. However, large-scale problems may quickly exceed the memory capacity even of very large machines, since the matrices are fully populated. We aim at illustrating how ℌ-matrices with adaptive cross approximation can provide an effective method to increase the size of the largest solvable problems by means of boundary element methods based on stream functions with modest implementation effort. The method is first validated on a benchmark problem for which an analytical solution is available and then applied to a complex problem of engineering interest related to the computation of the currents induced in the stabilizing shell of a nuclear fusion device.