MATLAB: Implementing Least Squares Estimator for a Given Model
The equation you're trying to solve is
$$
\mathbf{y}=\mathbf{X}\mathbf{h},
$$
where $\mathbf{h}$ is your unknown. The matrix $\mathbf{X}$ is going to have a time-shifted structure that reflects the convolution operator. If we assume that the $\mathbf{y}$ vector starts with y(3) i.e. ignores the first two zeroed out elements of y,  then the corresponding $\mathbf{X}$ matrix is given by:

$$
\mathbf{X}=\left[\begin{matrix}
x(3) &amp; x(2) &amp; x(1)\\
x(4) &amp; x(3) &amp; x(2) \\
x(5) &amp; x(4) &amp; x(3) \\
... &amp; ... &amp; ...
\end{matrix}  \right]
$$

You can add as many rows as you are have observations to support them. You can then solve for $\mathbf{h}$ by 
$$
\mathbf{h}_{est}=\text{pinv}(\mathbf{X})\mathbf{y}
$$