Plotting quiver plots in matplotlib

I want to plot the slope field for: 0.5*sin(0.5*pi*x)*sqrt(y+7)

import numpy as np
import matplotlib.pyplot as plt

# Specify the grid of dots 
x = np.arange(-3,3,0.3) 
y = np.arange(-2,4,0.3)
X, Y = np.meshgrid(x,y)

# Create unit vectors at each dot with correct slope
dy = 0.5*(np.sin(x*np.pi*0.5))*np.sqrt(y+7)
dx = np.ones(dy.shape)
norm = np.sqrt(X**2 + Y**2) 
dyu = dy/norm
dxu = dx/norm

# Plot everything 
plt.quiver(X,Y,dxu,dyu,)
plt.show()

and I got the second image below. I am trying to replicate the first image. How can I make uniform slope lines like that? Why am I getting variable-length lines anyway?

So, the first thing to do is remove the arrowheads, which @jasonharper shows in the comments can be done by adding these options to your quiver call: headwidth=0, headlength=0, headaxislength=0.

Next is to deal with the length. You're currently normalizing by X and Y when you should be normalizing by dx and dy. I would actually redefine dx and dy to use X and Y so they have the same shapes. And you will also want to make the pivot="mid" change that @jasonharper mentioned in their comment. With these changes, your code would look like this:

import numpy as np
import matplotlib.pyplot as plt

plt.close("all")

x = np.arange(-3, 3, 0.3) 
y = np.arange(-2, 4, 0.3)
X, Y = np.meshgrid(x, y)

dX = np.ones_like(Y)
dY = 0.5*np.sin(X*np.pi*0.5)*np.sqrt(Y+7)
norm = np.sqrt(dX**2 + dY**2)
dXu = dX/norm
dYu = dY/norm

plt.quiver(X, Y, dXu, dYu,
           headwidth=0, headlength=0, headaxislength=0, 
           pivot="mid")
plt.show()

Result:

To control the length, you'll need to change two things. The first thing would be to multiply dX and dY by the desired length. The second thing would be to add the following options to quiver: angles="xy", scale_units="xy", scale=1 to make sure the vectors are the desired length.

import numpy as np
import matplotlib.pyplot as plt

plt.close("all")

x = np.arange(-3, 3, 0.3) 
y = np.arange(-2, 4, 0.3)
X, Y = np.meshgrid(x, y)

dX = np.ones_like(Y)
dY = 0.5*np.sin(X*np.pi*0.5)*np.sqrt(Y+7)
norm = np.sqrt(dX**2 + dY**2)
dXu = dX/norm
dYu = dY/norm

L= 0.2
dXs = L*dXu
dYs = L*dYu

plt.quiver(X, Y, dXs, dYs,
           headwidth=0, headlength=0, headaxislength=0, 
           angles="xy", scale_units="xy", scale=1,
           pivot="mid")
plt.show()