import matplotlib.pyplot as plt
import numpy as np
import math as m
N=1000
a=1
b=2
nbins=25

rn=np.zeros(N)
rndm=np.random.rand(N)
i=0
for x in rndm:
    n=a+m.sqrt((2*a*rndm[i])/b)
    rn[i]=n
    i=i+1
xmin=rn.min()
xmax=rn.max()
fig,ax=plt.subplots()
plt.axvline(x=b,ymax=(N*(b-a)/nbins*(b/a*b-2))/100,c='red')
x=np.linspace(a,b)
y=N*(b-a)/nbins*(b/a*x-2)
ax.plot(x,y,c='red')
x=np.linspace(0,a)
nbins=np.linspace(a,b,nbins)
ax.hist(rn,bins=nbins,histtype='step',range=[0,b+1])
plt.xlim(0,b+1)
plt.ylim(0,100)
plt.show()