Firstly, I will expand on my comment.
If $a\neq 0$, then $0=-ay+e^{-y}$, which can be rearranged to $ye^y=\dfrac{1}{a}$. This has the "analytic" solution using the LambertW function $y=LambertW(1/a)$.
If $a=0$, then the original ODE is $y'=e^{-y}$, which has solution $y(x)=\log(x+c)$. This does not have a steady state.
Now, given that you are in a numerical analysis class, you can also solve this equation using numerical methods. Hence, a solution could be obtained using Newton's method. Given $f(y)=e^{-y}-ay$, the iteration would be
$$
y_{n+1}=y_n-\frac{e^{-y_n}-ay_n}{-e^{-y_n}-a}=y_n+\dfrac{1+ay_ne^{y_n}}{1+ae^{y_n}}.
$$
This will break down when $1+ae^y=0$.
Hence, if $a>0$, The iteration will not break down and is guaranteed to find a solution.
If $a<0$, it is possible that the expression is satisfied. Now, appealing to the geometric interpretation of the problem, we see that there is a solution provided that the graphs of $l_1=ay$ and $l_2=e^{-y}$ intersect. The least magnitude negative value for $a$ must then be the value where these two graphs are tangent. This leads to $a=-e^{-p}$ where $y=p$ is the point of tangency. But, we also require the line $l_1=ay$ to be satisfied at $(p,l_2(p))$ Hence, $e^{-p}=-e^{-p}p\implies p=-1$. Thus, the least magnitude negative value for $a$ is $a=-e$. The method will then also not breakdown if $a\leq-e$. This result is confirmed as the $LambertW(x)$ function is defined only when $x\geq-1/e$, which corresponds to $a\leq-e$ when $a<0$.
Conclusion
A solution will exist provided $a>0$ or $a\leq-e$. This is confirmed by this implicit plot of the solution.
The homogeneous form of the solution is actually
$$X_H=c_1e^{-t}sin(5t)+c_2e^{-t}cos(5t)$$
which exponentially decays, so the homogeneous solution is a transient. The steady state solution is the particular solution, which does not decay.
Best Answer
In different areas, steady state has slightly different meanings, so please be aware of that.
We want a theory to study the qualitative properties of solutions of differential equations, without solving the equations explicitly.
Moreover, we often want to know whether a certain property of these solutions remains unchanged if the system is subjected to various changes (often called perturbations).
It is very important to be able to study how sensitive the particular model is to small perturbations or changes of initial conditions and of various paramters.
This leads us to an area of DEQ called Stability Analysis using phase space methods and we would consider this for both autonomous and nonautonomous systems under the umbrella of the term equilibrium.
Autonomous
The equilibrium solution ${y_0}$ is said to be unstable if it is not stable.
Equivalent definitions can be written for the nonautonomous system $y' = f(t, y)$.
Now we can add notions of globally asymptoctically stable, regions of asymptotic stability and so forth.
From all of these definitions, we can write nice theorems about Linear and Almost Linear system by looking at eigenvalues and we can add notions of conditional stability.
Update
You might also want to peruse the web for notes that deal with the above. For example DEQ.
Regards