Using a Projection matrix to project grid

Question

I'm trying to project a grid using a projection matrix creating a grid. I generate some lines given two points. Then, I found a projection matrix and find the point images but when I'm trying to plot the result, I obtain weird things

def plot_grid(grid_range = [-10., 10.], num_lines = 10, return_points = False):
    # generate points
    lines_x = []
    lines_y = []
    factor = sum(np.abs(grid_range)) / num_lines

    for i in range(num_lines + 1) :
        
        init  = [grid_range[0] + factor * i, grid_range[0], 0., 1.]
        final = [grid_range[0] + factor * i, grid_range[1], 0., 1.]
        lines_x.append([init, final])

        init  = [grid_range[0] , grid_range[0] + factor * i, 0., 1.]
        final = [grid_range[1], grid_range[0] + factor * i, 0., 1.]

        lines_y.append([init, final])

    lines_x = np.array(lines_x)
    lines_y = np.array(lines_y)

    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d import Axes3D
    
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    for pair in lines_x :
        ax.plot(pair[:, 0], pair[:, 1], pair[:, 2], '-b')

    for pair in lines_y :
        ax.plot(pair[:, 0], pair[:, 1], pair[:, 2], '-b')
    ax.set_xlabel("x label")
    ax.set_ylabel("y label")
    ax.set_zlabel("z label")

    plt.show()
    
    
    if return_points :
        return lines_x, lines_y
    else :
        return


# Projection Matrix
image_width  = 640
image_height = 700
# Example of perspective projection
C      = np.zeros((3,1))
C[0], C[1], C[2] = 0, 0.5, 1
K      = np.eye(3)
focal  = 450
K[0,0] = focal
K[0,2] = image_width
K[1,1] = focal
K[1,2] = image_height
K[2,2] = 1.

theta   = 0
R       = np.zeros((3,3))
R[0][0] = np.cos(theta)
R[1][1] = np.cos(theta)
R[0][1] = -np.sin(theta)
R[1][0] = np.sin(theta)
R[2][2] = 1


P = K@R@np.concatenate([np.eye(3),-C],axis=1)
print(P)


# projected points
plt.axis('equal')
plt.xlim((0,image_width))
plt.ylim((0,image_height))
plt.gca().invert_yaxis()

for pair in lines_x :
    # Define a 3D square
    # Form its image
    Images = P@pair.T
    Images = Images.T
    # homogeneous coordinates
    p1 = [Images[0][0] / Images[0][2], Images[0][1] / Images[0][2]]
    p2 = [Images[1][0] / Images[1][2], Images[1][1] / Images[1][2]]
    plt.plot(p1, p2, '-b')
    
for pair in lines_y :
    # Define a 3D square
    # Form its image
    Images = P@pair.T
    Images = Images.T
    # homogeneous coordinates
    p1 = [Images[0][0] / Images[0][2], Images[0][1] / Images[0][2]]
    p2 = [Images[1][0] / Images[1][2], Images[1][1] / Images[1][2]]
    plt.plot(p1, p2, '-b')
    
plt.show()

but I obtain weird this like this

are something wrong with the parameters I'm choosing?