Here is an example of application btw. (Maybe others can find it useful). It works but I think it’s doing one full redundant backward pass.
def dfdx(x,param,f):
    x = Variable(x.data,requires_grad=True)
    z = f(x,param)
    z.backward([torch.ones_like(z)],create_graph=True)
    dfdx = x.grad

    return dfdx

class Model(torch.nn.Module):
    def __init__(self):
        super(Model, self).__init__()
        self.linear1 = nn.Linear(1, 1)
        self.linear2 = nn.Linear(1, 1)
        self.functional = lambda x,y : 100+(x+y)**2

    def forward(self, x):
        x = self.linear1(x)
        z = self.linear2(x).tanh() # Should be learned to be 0
        # this functional should be 2(x+y) if diff works. 
        x = dfdx(z,x,self.functional) 
        # x = dfdx(x,z,self.functional) # Both works. 
        return x

model = Model()

criterion = nn.MSELoss()
optimizer = torch.optim.Adam(model.parameters(), lr=1e-3)

n = 5
losses = []
x = Variable(torch.ones(n,1).cumsum(0)) #1,2,3,4,5
for i in range(1000):
    def closure():
        x = Variable(torch.randn(n,1)) 
        y = model(x)
        loss = criterion(y,x)
        optimizer.zero_grad()
        loss.backward()
        print(loss.data.numpy())
        return loss    
    optimizer.step(closure)
print(model(x))