Relationship between Discrete Deconvolution and Toeplitz Matrices
I think you are constructing the convolution matrix incorrectly. Let's do an easy example with some numbers to see how it is constructed

Let's imagine that youhave the following vectors of dimensions $M$ and $N$, And let's for the sake of it, say that $M = 5$ and $N = 4$

$\mathbf{a} = [a_{1},a_{2},a_{3},a_{4}]^{T}$

$\mathbf{b} = [b_{1},b_{2},b_{3},b_{4},b_{5}]^{T}$

$
\mathbf{c} = \mathbf{a}*\mathbf{b}
$

In matrix form, if we assume $\mathbf{a}$ is the system response, the convolution is given by

\begin{equation} 
    c = \mathbf{A}\mathbf{b}
\end{equation}

with

\begin{equation} 
    \mathbf{A} = \left[
    \begin{array}{ccccc}
        a_{1} &amp; 0 &amp; 0 &amp; 0 &amp; 0 \\
        a_{2} &amp; a_{1} &amp; 0 &amp;  0 &amp; 0 \\
        a_{3} &amp; a_{2} &amp; a_{1} &amp;  0 &amp; 0 \\
        a_{4} &amp; a_{3} &amp; a_{2} &amp; a_{1} &amp;  0 \\
        0 &amp; a_{4} &amp; a_{3} &amp; a_{2} &amp; a_{1} \\
        0 &amp; 0 &amp; a_{4} &amp; a_{3} &amp; a_{2} \\
        0 &amp; 0 &amp; 0 &amp; a_{4} &amp; a_{3} \\
        0 &amp; 0 &amp; 0 &amp; 0 &amp; a_{4} \\
    \end{array}\right]
\end{equation}

By doing the convolution with this matrix, you get the full operation including the transients. So the correct matlab code would be

A = toeplitz([a,zeros(1,M-1)],[a(1),zeros(1,M-1)]) 

(Check help toepliz to verify how it constructs the matrix)

Use this matrix and the result of the convolution is exactly as with the matrix multiplication.

Now, a word of wisdom, as you can clearly see, the full convolution matrix with the transients is not square, so in order to "invert" it and do deconvolution you will need to use the Moore-Penrose Pseudoinverse. And you will get some differences in the beginning and end of the resulting vector, but it is the best you can do, as using the pseudoinverse would be the same as performing least squares, since the deconvolution is given by

$\mathbf{b} = \left(\mathbf{A}^{T}\mathbf{A}\right)^{-1}\mathbf{A}^{T}\mathbf{c}$

This result will always hold, because the dimensions of $A$ are always $(M+N-1)\times M$, so the matrix is always "skinny" and the equation holds.

Try this, it should give you the result you were looking for.