Induction hardening is a well known method for the heat treatment of steel components. With the concept of multi-frequency hardening, where currents with two different frequency components are provided on a single inductor coil, it is possible to optimize the hardening zone to follow a given contour, e.g. of a gear. In this article, we consider the simulation of multi-frequency induction hardening in 3D. The equations to solve are the vector potential formulation of Maxwell's equations describing the electromagnetic fields, the balance of momentum to determine internal stresses and deformations arising from thermoelasticity and transformation induced plasticity, a rate law to determine the distribution of different phases and the heat equation to determine the temperature distribution in the workpiece. The equations are solved using adaptive finite element methods. The simulation results are compared to experiments for discs and for gears. A very good agreement for the hardening profile and the temperature is observed. It is also possible to predict the distribution of residual stresses after the heat treatment