The finite difference element method (FDEM) is a black-box solver for the solution of nonlinear systems of elliptic and parabolic PDEs. An algorithm has been developed to generate on an unstructured FEM grid difference formulas of arbitrary consistency order q. From the difference of difference formulas of different consistency order, an estimate of the discretization error is obtained. An error equation permits the explicit following of all errors and gives the prescriptions for the selfadaptation of the method. Coupled domains with different PDEs and different nonmatching grids that slide relative to each other can be treated and a global error estimate is computed. Thus, we get an FDM that is in all aspects more flexible than the FEM. The whole code is efficiently parallelized on distributed memory parallel computers.