What I'm trying to do is simply a variation of a heavyside function, for 3 before the discontinuity and 4 after. But when I plot it, at the discontinuity it's not a straight vertical line but a sort of $x=\alpha*x$ function. Can someone explain why it's like that?
What happens as you add more points to your domain near $x=x_c$? Also, many (most, all(?)) definitions of the Heaviside step function give a value of the mean of the extremes at the transition point, but you initialize to zero and skip $x=x_c$ in your if-block. You should handle that case, too.
Also, why write your own when numpy supplies one?
Edited to address the comment about the Python numpy.heaviside() function (you can have exactly what you ask for):
import numpy v=numpy.linspace(-1,1,num=21) out=2.*numpy.heaviside(v,0.5)+3. array([3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 4., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5.])
Here I've carefully chosen my parameters to get a point right on the transition point so that there's one there on the output. As for any value you want, you can always manipulate the result to adjust it to the domain and range you like.