A variational time discretization of anisotropic Willmore flow combined with a spatial discretization via piecewise affine finite elements is presented. Here, both the energy and the metric underlying the gradient flow are anisotropic, which in particular ensures that Wulff shapes are invariant up to scaling under the gradient flow. In each time step of the gradient flow a nested optimization problem has to be solved. Thereby, an outer variational problem reflects the time discretization of the actual Willmore flow and involves an approximate anisotropic L 2 -distance between two consecutive time steps and a fully implicit approximation of the anisotropic Willmore energy. The anisotropic mean curvature needed to evaluate the energy integrand is replaced by the time discrete, approximate speed from an inner, fully implicit variational scheme for anisotropic mean curvature motion. To solve the nested optimization problem a Newton method for the associated Lagrangian is applied. Computational results for the evolution of curves underline the robustness of the new scheme, in particular with respect to large time steps.