Hi Alex, thanks for sharing this link. It helped me finding out where I did go wrong and at the end I got the same result with my implementation as with the jacobian function from pytorch.
So in case anybody is reading this later on: My fault was that I calculated the gradient with respect to the normalized input but I should have done it with respect to the raw unnormalized input. Therefore I needed to add another function which normalizes the raw input and keeps the gradient. So here is the correct implementation:
def calculate_jacobi(model, inp_data, inp_scaler, out_scaler):
	# set requires_grad to True so that the gradient is calculated with respect to unnormalized input
	inp_data.requires_grad = True
	
	normalized_input = normalize_inp_keeping_gradient(inp_data, inp_scaler)

	# set requires_grad for all network parameters to False so that the gradient is not
	# calculated w.r.t. the network parameters
	for param in model.parameters():
	    param.requires_grad = False

	output = model(normalized_input)
	unnorm_output = revert_output_normalization_keeping_gradient(output)

	# calculate gradient for every output dimension -> jacobi matrix
	jacobi_matrices = []
	for out in unnorm_out:
		jacobi_matrix = []
		for out_dim in out:
			out_dim.backward(retain_graph=True)
			# calculate gradient with respect to unnormalized input
			jacobi_matrix.append(inp_data.grad.numpy())
			inp_data.grad = None
		jacobi_matrices.append(np.array(jacobi_matrix))

	return np.array(jacobi_matrices)

def normalize_inp_keeping_gradient(inp, scaler):
	normalized_inp = torch.zeros(inp.shape)
	for dim in range(0, len(normalized_inp[0])):
		minimum = scaler.min_[dim]
		scale = scaler.scale_[dim]
		normalized_inp[:, dim] = torch.add(torch.mul(inp[:, dim], scale), minimum)

	return normalized_inp

def revert_output_normalization_keeping_gradient(normalized_output, scaler):
	unnorm_output = torch.zeros(normalized_output.shape)
	for dim in range(0, len(unnorm_output[0])):
		minimum = scaler.min_[dim]
		scale = scaler.scale_[dim]
		unnorm_output[:, dim] = torch.div(torch.sub(normalized_output[:,dim], minimum), scale)

	return unnorm_output

And corresponding to the implementation above one could use the jacobian function of pytorch as follows:
example_input = torch.FloatTensor([[5.543]])
torch.autograd.functional.jacobian(func_to_derive, example_input, create_graph=True))

where the func_to_derive is defined as follows:
def func_to_derive(model, inp_data, inp_scaler, out_scaler):
	inp_data.requires_grad = True

	# normalize input with keeping its gradient
	normalized_inp_data = inp_scaler.normalize_inp_keeping_gradient(inp_data)

	# set requires_grad for all network parameters to False so that the gradient is not
	# calculated w.r.t. the network parameters
	for param in model.parameters():
		param.requires_grad = False

	# calculate corresponding output
	normalized_output = model(normalized_inp_data)

	# revert normalization for output so gradient of unnormalized output can be calculated
	unnorm_out = out_scaler.revert_output_normalization_keeping_gradient(normalized_output)

	return unnorm_out