Your ray direction calculation is quite odd. Firstly, it's not normalized which will cause you to overshoot/undershoot the surface in most cases. Secondly, the projection is very odd - not sure what it is but it's definitely not correct. Maybe it was supposed to be the direction from (0, 0, 0) to the calculated point?

Here's some code I wrote for a very similar raytracer. There are definitely smarter and more customizable (rotation, FOV, etc.) ways of calculating this but I came up with this and it gives me good results, at least for a quick test render. It's in python but you should still be able to get the gist of it. It assumes that the camera is facing along the -Z axis.

# x, y, width, and height are in pixels
u = x / width 
v = y / height

top_left     = Vec3(-1, +1, -1)
top_right    = Vec3(+1, +1, -1)
bottom_left  = Vec3(-1, -1, -1)
bottom_right = Vec3(+1, -1, -1)
# Note `1 - v` to flip to the correct way up
direction = vec.bilerp(
    top_left, top_right, bottom_left, bottom_right, u, 1 - v
)

where vec.bilerp is something like

    lerp(
        lerp(top_left, top_right, u),
        lerp(bottom_left, bottom_right, u),
    v)

You might not need to flip it with 1 - v depending on how you output your image. Also, I just noticed that my code doesn't normalize the direction either. *facepalm*