Hey,
I think the confusion comes from what the backward pass computes.
if you function is y = relu(x) and you have a loss L.
Then the backward is computing dL/dx = dL/dy * dy/dx. Where dy/dx is the jacobian of the relu function. But keep in mind that you need to do a matrix product with the gradient flowing back.
The Jacobian of the relu is very close to an identity matrix with 1 on the diagonal for positive inputs and 0 on the diagonal for negative inputs.
So the matrix product here will effectively 0 out the entries in the grad flowing back when the input was negative. And just keep it as is when the input was positive.
In the tutorial code,
    grad_h = grad_h_relu.copy()

Is here only to make sure you don’t modify inplace the gradient flowing back inplace (in case it would be used somewhere else.
    grad_h[h < 0] = 0

Zero out the entries corresponding to inputs that were negative.
Note that this is single op is the same as doing the matrix product from the chain rule.
In your code sample, grad = x.copy() does not look right. x should be input to the forward pass while grad should be the gradient flowing back (the input of the backward function).