What does the descision boundary of a relu look like?

What would be the shape of this separator (decision boundary) in case we also take a relu on the output and only then threshold?

For just a single neuron, indeed the decision boundary will just be a plane.

I keep thinking it has to be non linear, otherwise NNs wouldn't work...

Single-layer neural networks don't demonstrate the same nonlinearity as multi-layer networks, only as much as an activation function. (With sigmoid activation, you're just doing logistic regression.)

As a follow up, I would like to understand multi level descision boundaries...

With multiple layers, you gain the ability to combine boundary planes into polytopal regions. A ReLU-only neural network represents a piecewise-linear function, with finitely many polytopal (i.e., flat-edged) pieces. That's clear for every neuron in the first layer. In each subsequent layer, you take a linear combination of such functions, which is again such a function; applying a ReLU to that then can cut each region in two, but the result is still such a function. (You may also be interested in this paper, and this expository blog post.)

You can see some of this in the tensorflow playground. In a version of that link, I've grabbed the weights and put the resulting function into Wolfram|Alpha. You can see how the boundary lines of each first-layer neuron shows up as a kink, and each of the resulting regions supports a linear function. Then slicing at a fixed height gives you the hexagonal region that works well enough to classify the circle dataset. Increasing the number of layers or neurons per layer allows you to cut the space into more regions.

Then, even if the output neuron applies a sigmoid as in your tensorflow playground, the decision boundary is just a level set of a piecewise-linear function with all linear boundaries, which is itself linear.