How to find center of a conic section from the equation?
Matrix derivation
The original equation can be described in a more concise form as $$\vec x A \vec x^T + \vec x \vec b^T + \vec b \vec x^T + f =0$$ or $$\vec x A \vec x^T + 2\vec x \vec b^T + f =0$$ where $\vec b=(d,e)$ and $A=\begin{pmatrix}a&b\\b&c\end{pmatrix}$.
If the vector $\vec x_0$ represents the center then in a coordinate system with the origin in this point the equation will be $$(\vec x-\vec x_0) A (\vec x-\vec x_0)^T + g =0$$ which can be rewritten as $$\vec x A \vec x^T - 2\vec x A\vec x_0^T + \vec x_0 A \vec x_0^T +g.$$
Which means the we need to choose $\vec x_0$ in such way that $-2A\vec x_0^T=2\vec b^T$, i.e. $$A\vec x_0^T=-\vec b^T,$$ $$\begin{pmatrix}a&b\\b&c\end{pmatrix} \begin{pmatrix}x\\y\end{pmatrix}= \begin{pmatrix}d\\e\end{pmatrix}.$$
This is precisely the following system of linear equations written in the matrix form: $$ \begin{align*} ax + by &= -d\\ bx + cy &= -e \end{align*} $$
(To give credit where credit is due, I was shown this derivation by a colleague of mine who taught labs for the same subject as I did, but with another group of students.)
Derivatives
We can view the curve as a level set of the function $$F(x,y)=ax^2+2bxy+cy^2+2dx+2ey.$$ The level sets of this function create the family of curves $ax^2+2bxy+cy^2+2dx+2ey+f=0$, where the parameter $f$ is changing.
For example, if $\delta>0$ we get concentric ellipses, like in this example - the plots are taken from WolframAlpha:
And a different plot:
If $\delta>0$, we get a family of hyperbolas. WolframAlpha:
And a different plot:
In the first case, the function $F$ has minimum at the center of each of the ellipse. In the second case there is a saddle point. In both cases, the center fulfills $$\frac{\partial F}{\partial x}= \frac{\partial F}{\partial y} =0$$ which leads to the system of equations \begin{align*} ax + by +d&= 0\\ bx + cy +e&= 0 \end{align*} which is precisely the linear system described above.
So far we have only drawn some pictures (which is good for intuition, but picture is not a proof.) If we want to give a more rigorous reason why this works we can use some known facts about quadratic forms. We know from Principal axis theorem that there is a rotation and translation such that in the new coordinates the curve has $$\lambda_1 x^2 + \lambda_2 x^2=const$$ Clearly $F(x,y)=\lambda_1 x^2 + \lambda_2 x^2$ has partial derivatives equal to zero at the origin. (And depending on the signs of $\lambda_{1,2}$ we can say whether it has minimum, maximum or saddle point there.)
I would emphasize that a pure translation in the plane preserves the relationship between a figure and its center. Given $$ A x^2 + 2 B xy + C y^2 + 2 D x + 2 E y + F = 0 $$ and solving for $(x_0, y_0)$ in $$ A x_0 + B y_0 + D = 0, $$ $$ B x_0 + C y_0 + E = 0, $$ we introduce translated coordinates $(u,v)$ with $$ x = u + x_0, $$ $$ y = v + y_0. $$ I'm not sure how much I will type in, but the original equation
$$ A x^2 + 2 B xy + C y^2 + 2 D x + 2 E y + F = 0 $$ becomes $$ A u^2 + 2 B uv + C v^2 = \rm{constant.} $$ Regardless of what type of figure this might be, it is centrally symmetric in $(u,v)$ coordinates: if $(u,v)$ is a solution, so is $(-u,-v).$
APPENDIX:
$$ A x^2 + 2 B xy + C y^2 + 2 D x + 2 E y + F = 0. $$
$$ x = u + x_0 $$ $$ y = v + y_0 $$
$$ x^2 = u^2 + 2 x_0 u + x_0^2 $$ $$ xy = uv + y_0 u + x_0 v + x_0 y_0 $$ $$ y^2 = v^2 + 2 y_0 v + y_0^2 $$
$$ A x^2 + 2 B xy + C y^2 = (A u^2 + 2B uv + C v^2) + (2Ax_0 + 2 B y_0) u + (2Bx_0 + 2 C y_0) v + (Ax_0^2 + 2B x_0 y_0 + C y_0^2) $$ $$ 2 D x + 2 E y + F = 2 D u + 2 E v + 2 D x_0 + 2 E y_0 + F $$ In the original $A x^2 + 2 B xy + C y^2 + 2 D x + 2 E y + F,$ the quadratic terms come out to $A u^2 + 2B uv + C v^2,$ as was guaranteed to happen for the highest degree terms. Next, we get $ (2Ax_0 + 2 B y_0 + 2 D)u, $ so that the coefficient of the linear term $u$ is actually $0.$ Also, $ (2Bx_0 + 2 C y_0 + 2 E)v , $ so that the coefficient of the linear term $u$ is actually $0.$ The constant term comes out to $$ Ax_0^2 + 2B x_0 y_0 + C y_0^2 + 2 D x_0 + 2 E y_0 + F $$ to which we say, Why Not? The whole thing is now $$ A u^2 + 2B uv + C v^2 + ( Ax_0^2 + 2B x_0 y_0 + C y_0^2 + 2 D x_0 + 2 E y_0 + F) = 0, $$ where the thing in parentheses is a constant. Let us give the constant a name, $ W = ( Ax_0^2 + 2B x_0 y_0 + C y_0^2 + 2 D x_0 + 2 E y_0 + F). $ So we get $$ A u^2 + 2B uv + C v^2 + W = 0 $$