Tracing the curve of maximum slope on a multivariable function given initial conditions. Let's say we have the multivariable function $f(x,y)=-x^2-xy-y^2$. Then $\nabla f=\frac{\partial f}{\partial x} \hat{i}+\frac{\partial f}{\partial y}\hat{j}=(-2x-y)\hat{i}+(-x-2y)\hat{j}$.
$\nabla f$ gives the direction vector where the function's slope is highest.
Here, I've plotted the function on Wolfram Mathematica:

What I want to do is as follows:

Given an object to follow the vector with maximum slope and move along the surface $f$, I want to obtain an equation for the curve traced by that object (which I've sketched in $\color{red}{\text{red}}$), given some initial point. Let's say we start at the point $(x,y)=(10,10)$.

Below is my full attempt to do this:
I've calculated the point at which the maximum value of $f$ is by letting each partial derivative equal to zero, and solving the system. I have done this to verify whether the traced line I find will stop at the origin $(0,0,0)$.
$$\frac{\partial f}{\partial x}=-2x-y=0$$
$$\frac{\partial f}{\partial y}=-x-2y=0$$
Which gives me the point $(0,0)$, and thus the maximum value is $f=0$.
To solve for the curved line, I've thought of parametrizing the curve with respect to some variable $t$.
Therefore, I've made the following system of ODE's (From $\nabla f$):
$$\frac{dx}{dt}=-2x-y$$
$$\frac{dy}{dt}=-x-2y$$
With initial conditions $x(0)=10$ and $y(0)=10$.
From this, I've solved it using matrices and finding eigenvalues and eigenvectors:
$$\mathbf{Y'}=\begin{pmatrix} -2 & -1 \\ -1 & -2 \end{pmatrix} \mathbf{Y}$$
Where $\mathbf{Y}=\begin{pmatrix} x \\ y \end{pmatrix}$.
Using $\det(A-\lambda I)=0$, I got:
$$\begin{vmatrix} -2-\lambda & -1 \\ -1 & -2-\lambda \end{vmatrix}=(-2-\lambda)^2-1=\lambda^2+4\lambda+3=0$$
Solving this quadratic gives me the eigenvalues $\lambda_1=-3$ and $\lambda_2=-1$.
The corresponding eigenvectors are thus:
$$\vec{v_1}=\begin{pmatrix} 1 \\ 1 \end{pmatrix} \qquad \vec{v_2}=\begin{pmatrix} -1 \\ 1 \end{pmatrix}$$
Therefore, I've obtained the following general solution:
$$\begin{pmatrix} x(t) \\ y(t) \end{pmatrix}=c_1 \begin{pmatrix} 1 \\ 1 \end{pmatrix}e^{-3t}+c_2 \begin{pmatrix} -1 \\ 1 \end{pmatrix}e^{-t}$$
When substituting initial conditions for $(x(0),y(0))=(10,10)$ I obtained the system:
$$\begin{cases} 10=c_1-c_2 \\ 10=c_1+c_2 \end{cases}$$
Which gave me $c_1=10$ and $c_2=0$.
Therefore, the specific solution for this case is:
$$x(t)=10e^{-3t}$$
$$y(t)=10e^{-3t}$$
I figured that I could substitute the solutions into $f(x,y)$ to obtain an expression for $f(t)$.
$$f(t)=-(10e^{-3t})^2-(10e^{-3t} \cdot 10e^{-3t})-(10e^{-3t})^2=-300e^{-6t}$$
Therefore, I figured that these parametrized curves would give me the red line. It seems correct to me, since as $t \to \infty$, $x=y=f=0$, which is what was predicted.
If I reverse the parametrization, I get the following curve (with restricted domain, of course):
$$x^2=y^2=-\frac{30000}{z}$$
However, I have several questions since I've made up this method myself (Even though it probably already exists):

1) Is my attempt correct? I am wondering since I am unsure whether the system of equations is actually correct to begin with.
2) Will the red line always be a geodesic (shortest route) along the surface given any initial conditions? If not, is it for the initial condition $(10,10)$?
3) Is there a simpler method to do this? Is there a more general equation to do this? I am aware of the geodesic equation (below), however I am unsure how to solve it, or if it is relevant at all.
$$\frac{d^2 \gamma^\lambda}{dt^2}+\Gamma^{\lambda}_{\mu \nu} \frac{d\gamma^{\mu}}{dt} \frac{d\gamma^{\nu}}{dt}=0$$

Please let me know if there are any problems with the questions. Any thoughts on the problem will be appreciated. Thanks in advance.
 A: I've double checked your solution and it is ok. The steps are ok and I have drawn the curve and it is exactly the expected path: you chose a starting point laying into one of the surface's planes of symmetry, so the curve have to follow that plane into the surface. I don't know a simpler way different from what you chose in the general case. Surely there is not a simpler way (I have seen these curves before: these "curves of maximum slope", better said, their projection on the xy-plane, are identificable as the "field lines" associated to the function f understood as some kind of potential)
I have done my own "experiments", by choosin as starting point (10,0) (far from any plane of symmetry). The solution of your ODE is now:
$$x(t)=5e^{-3t}+5e^{-t}$$
$$y(t)=5e^{-3t}-5e^{-t}$$
$$z(t)=-25\left(3e^{-6t}+e^{-2t}\right)$$
This curve reaches the maximum of f too (all curves do) with some surprise by me and with no obvious clue about whether or not it's a geodesic. Further, I don't find another way to check in the general case wich one of these "maximum slope curves" are geodesics other than checking your general solution into the geodesic equations you have quoted (it takes a little amount of patience calculate all those $\Gamma's$) It is interesting for me to know that. Of course, it is true that the family of geodesics is "larger" than the one of "maximum slope"
