The key relation for functional differentiation is $$\frac{\delta}{\delta f(y)}f(x)=\delta(x-y), $$ where $\delta(x-y)$ is the Dirac delta function, and the usual properties of differentiation (e.g. linearity, chain rule) still hold. Another example: $$\frac{\delta}{\delta f(y)}(f(x)f(x))=2 \delta(x-y) f(x).$$ A possible Mathematica version could look like this FunctionalD[f[x],f[y]]=DiracDelta[x-y] How would one go about implementing this?