Deriving the average depth for a randomly generated binary search tree
Suppose that $T$ is a random tree on $n$ vertices, according to the distribution you consider. Let $T_L,T_R$ be its two subtrees, and let their size be $N_L,N_R$. Let $P$ be the total path length of $T$. Let $P_L(i)$ be the total path length of $T_L$ if $N_L=i$, and zero otherwise, and define $P_R(i)$ similarly. We have
$$
P = \sum_{i=0}^{n-1} (P_L(i) + P_R(n-1-i)) + n-1.
$$
According to linearity expectation, we have 
$$
\mathbb{E}[P] = \sum_{i=0}^{n-1} (\mathbb{E}[P_L(i)] + \mathbb{E}[P_R(n-1-i)]) + n-1.
$$
By definition, $\mathbb{E}[P] = D(n)$, where $D(n)$ is the expected internal path length of a random tree on $n$ vertices.

According to your probability distribution, the probability that $N_L=i$ is $1/n$. When that happens, $T_L$ is just a random tree on $i$ vertices. Therefore
$$
\mathbb{E}[P_L(i)] = \Pr[N_L=i] \mathbb{E}[P_L(i)|N_L=i] + \Pr[N_L\neq i] \mathbb{E}[P_L(i)|N_L\neq i] = \\ \frac{1}{n} \cdot D(i) + \left(1-\frac{1}{n}\right) \cdot 0 = \frac{D(i)}{n}.
$$
Therefore we obtain the recurrence
$$
D(n) = \frac{1}{n} \sum_{i=0}^{n-1} (D(i) + D(n-1-i)) + n-1 = \frac{2}{n} \sum_{i=0}^{n-1} D(i) + n-1.
$$
Equivalently,
$$
nD(n) = 2\sum_{i=0}^{n-1} D(i) + n(n-1).
$$
The initial condition is $D(0) = 0$.

We can solve this recurrence using generating functions. Let $P(x) = \sum_{n=0}^\infty D(n) x^n$. First,
$$
P'(n) = \sum_{n=1}^\infty nD(n) x^{n-1}.
$$
Second,
$$
\frac{P(n)}{1-x} = \sum_{n=0}^\infty [D(0) + \cdots + D(n)] x^n.
$$
Third,
$$
\sum_{n=0}^\infty n(n-1) x^{n-2} = \frac{d^2}{dx^2} \frac{1}{1-x} = \frac{2}{(1-x)^3}.
$$
Putting everything together, we obtain
$$
P'(x) = \sum_{n=1}^\infty nD(n) x^{n-1} =
2\sum_{n=1}^\infty (D(0) + \cdots + D(n-1)) x^{n-1} +
\sum_{n=1}^\infty n(n-1) x^{n-1} = \\
\frac{2P(x)}{1-x} + \frac{2x}{(1-x)^3}.
$$
The solution to this ODE is
$$
P(x) = \frac{-2\log (1-x)-2x}{(1-x)^2} = \frac{1}{(1-x)^2} \sum_{n=2}^\infty \frac{2}{n} x^n,
$$
which implies that
$$
D(n) = 2\sum_{m=2}^n \frac{n-m+1}{m} = 2(n+1)(H_n-1)-2(n-1) = 2(n+1)H_n-4n.
$$
Since $H_n = \ln n + \gamma + 1/2n + O(1/n^2)$, we deduce that
$$
D(n) = 2(n+1)(\ln n + \gamma + 1/2n + O(1/n^2)) - 4n = \\
2n\ln n + 2\ln n + (2\gamma-4) n + 2\gamma + 1 + O(1/n).
$$