This is absolutely amazing! I never thought to use finite difference methods inside the loss function to compute derivatives. I tried to build PINNs in WL in the past, but I gave up after realizing that WL's neural network library based on mxnet could not automatically compute higher-order derivatives like TensorFlow or PyTorch. Your approach appears to be a promising path forward.
