How do I calculate the normals of a tubular mesh around any explicit function?

Question

I used this to help create a 3d cylindrical mesh that wraps around any y = f(x).
Each vertex is calculated by a point P:

  const P = (x,y,dy,r,u) => {
      const normal = -1/dy;
      return {
        x: x + r*Math.cos(u),
        y: y + r*Math.cos(u) * normal,
        z: r*Math.sin(u)
      }
    }

Images: x^2 sin(x)

As you can see there are issues when |dy| < 1. Obviously this is because of my calculation of the normal.
In general, I want the cylinders to face the proper direction, without any scaling.

I had tried to do a piecewise for |dy| < 1 but it isn't the elegant solution im looking for. I know that when the cylinder is facing vertically, the vertices are a function of x and z ( dy = inf so normal = 0 accounts for this), and when the cylinder is facing horizontally the vertices should be strictly in terms of y and z. This is where I get the divide by zero discontinuity and and the weird scaling when dividing by |dy| < 1

Basically, I am looking for some f(dy) that determines how much the x and y component contribute to the generation of the vertices.
full code