How does the evaluation of $\dot q_i$ at $q_1+\mathrm dq_1$ yield $\dot q_i +\dfrac{\partial \dot q_i}{\partial q_1}\mathrm dq_1 \;?$ I've been following Reif's Fundamentals of Statistical and Thermal Physics; there I came before the derivation of Liouville's theorem:

There I couldn't understood few things.
I could conceive the change in the number of systems in $\mathrm dt$ is given by $$\frac{\partial \rho}{\partial t}\; \mathrm dt\; (\mathrm d q_1,\mathrm dq_2, \ldots,\mathrm d q_f; \mathrm d p_1, \mathrm dp_2,\ldots, \mathrm d p_f)$$ where 


*

*$$\rho(q_1,q_2,\ldots, q_f; p_1,p_2,\ldots, p_f ; t)\;\mathrm d q_1,\mathrm dq_2, \ldots,\mathrm d q_f; \mathrm d p_1, \mathrm dp_2,\ldots, \mathrm d p_f = \textrm{no of systems in the ensemble at $t$ in the phase-space volume}\;(\mathrm d q_1,\mathrm dq_2, \ldots,\mathrm d q_f; \mathrm d p_1, \mathrm dp_2,\ldots, \mathrm d p_f)$$


But then, I couldn't understand why the number of systems 'entering this volume in time $dt$ through the face $q_1$= constant' is given by the quantity $\rho(\dot{q_1}\mathrm dt, \mathrm dq_2, \ldots,\mathrm dp_f )\;.$ 
My questions are:
$\bullet$ How does $\rho(q_1,q_2,\ldots, q_f; p_1,p_2,\ldots, p_f ; t)((\dot{q_1}\mathrm dt)(\mathrm dq_2, \ldots,\mathrm dp_f ))$ represent the number of systems that would enter the volume in time-interval $\mathrm d t\;?$
$\bullet$ How does the evaluation of $\dot q_i$ at $q_1+\mathrm dq_1$ yield $\dot q_i +\dfrac{\partial \dot q_i}{\partial q_1}\mathrm dq_1 \;?$
 A: It is simpler to think of the motions of various systems in the ensemble like that of the particles of a gas in a $2f$ dimensional space. Carrying the analogy further, let us try to answer the question how many gas particles moving with velocity $v_z$ would collide with a unit area of a wall perpendicular to say $z$-axis in unit time? The answer is, as many gas particles will collide as are contained in the volume $v_z$ (Remember unit area) and are moving with velocity $v_z$. If $\rho_{v_z}$ is the density of gas particles which are moving with velocity $v_z$, then the answer to the problem posed is evident $v_z\rho_{v_z}$. In the case of a Statistical Ensemble, the role of the velocities is played by the rates of change of coordinates of the phase space, i.e by $\dot{q_i}$ and $\dot{p_i}$ and gas density is replaced by the statistical distribution function.
As to your second query, note that $\dot{q_i}$ varies over the volume of the phase space, so the second term in the expression $\dot q_i +\dfrac{\partial \dot q_i}{\partial q_1}\mathrm dq_1 $ is to take account for such a "spatial" variation.
