Effect of Adding the cyclic prefix on the toeplitz matrix in OFDM
What you mean might be circulant matrix instead of toeplitz matrix. See section 3.4.4 in https://web.stanford.edu/~dntse/Chapters_PDF/Fundamentals_Wireless_Communication_chapter3.pdf about how the circular convolution in OFDM is represented by matrix operations (eq 3.130 onwards).

First, in almost all standard OFDM systems, you can assume $D \le L$. The cyclic prefix will be less than or equal to maximum multi-path delay, so as long as $D \le L$, the linear convolution $x * h$ gets converted circular convolution of $x$ and $h$. 

When $y_c = H x_{cp}(k)$, you only need to take $N$ original elements for $x$.  $H$ is size $N \times N$. This is because you have already re-written the circular convolution in matrix form. Each row of $H$ will do dot-product with $x$ to generate $y_c[n]$. Like that there will be $N$ values of $y_c$ corresponding to each row of $H$. 

$y$ is of length $N+L +D-1$. As you correctly mentioned $y_c = y(D+1:D+N)$.

Equivalent $H$ (size $N \times N$):

h(0)  0    0 ... h(L-1)  h(L-2) .. h(1)

h(1)  h(0) 0     0  ... h(L-1) h(L-2).. h(2)

h(2)  h(1)  h(0)  0 ..  0 0 h(L-1) h(L-2) .. h(3)

.

.

.

0 0 0 ...    h(L-1) h(L-2) .. h(2) h(1) h(0)

MATLAB command to generate above $H$:

Assuming $h$ is the vector representing channel having $L$ non-zero values.

H = toeplitz([h zeros(1,N-L)][h(1) zeros(1,N-L) flip(h(2:L))])


For a small example where $N=4, L=2, D=1$

$x = [x0, \,x1,\, x2,\, x3]$

$h = [h0\,,h1]$

$y = [x3h0\;, x3h1+x0h0\;, x0h1+x1h0\;, x1h1+x2h0\;, x2h1+x3h0\;, x3h1]$

$H$

h0 0 0 h1

h1 h0 0 0

0 h1 h0 0

0 0 h1 h0