In matrix terms the normal regression model can be written as
$$ Y=X\beta+\epsilon $$
where $E(\epsilon)=0$ and $\sigma^2\{\epsilon\}=\sigma^2I$, i.e., $\epsilon\sim N(0,\sigma^2I)$
or
$$ \begin{bmatrix} y_1\\y_2\\\vdots\\y_n \end{bmatrix}= \begin{bmatrix} 1 & x_1\\ 1 & x_2\\ \vdots & \vdots\\ 1 & x_n \end{bmatrix} \begin{bmatrix} \beta_0\\ \beta_1 \end{bmatrix} + \begin{bmatrix} \epsilon_1\\ \epsilon_2\\ \vdots\\ \epsilon_n \end{bmatrix} $$
where $\epsilon\sim N(0,\sigma^2)$
In short, find
$$ \beta=(X'X)^{-1}X'Y $$
We need to find $\beta$ minimize the quantity ($Q$, i.e., SSE or RSS):
$$ \begin{align*} Q&=\sum[y_i-(\beta_0+\beta_1 x_i)]^2\\ &=(Y-X\beta)'(Y-X\beta)\\ &=Y'Y-\beta'X'Y-Y'X\beta+\beta'X'X\beta \end{align*} $$
Since
Thus, we obtain
$$ Q=Y'Y-2\beta'X'Y+\beta'X'X\beta $$
Further,
$$ \begin{align*} \frac{\partial Q}{\partial \beta}&= \begin{bmatrix} \frac{\partial Q}{\partial \beta_0}\\ \frac{\partial Q}{\partial \beta_1} \end{bmatrix}\\ &=-2X'Y+2X'X\beta\\ -2X'Y+2X'X\beta&=-2X'Y+2X'X\beta\\ \beta&=(X'X)^{-1}X'Y \end{align*} $$