Solved by googlebot in post #2 
              
              
                You’d calculate & store dBessel(…)/dorder (generally, dF(x)/dx). This can be done in [autograd.Function’s] forward() or postponed to backward(). Irrespective of how you obtained a gradient, in backward you do multiplication: dLoss/dx = dLoss/dF(x) * dF(x)/dx = grad_out * local_grad  and return that …