What are the steps involved in solving a quartic polynomial modulo a prime modulus? This:
$$x^4 + 21x^3 + 5x^2 + 7x + 1 \equiv 0 \mod 23$$
Leads to:
$$x = 18 || x =19$$
I know this because of this WolframAlpha example and because a fellow member posted it in a since deleted & related question. 

What I don't understand are the steps involved in arriving at x = 18 || x = 19 from this equation.
My question starts with the reduced terms mod 23 example in the linked question. I'm now trying understand how to reduce this equation to x = 18 || x = 19.
I have come across a few posts and theorems that hint a solutions, but I lack the math skills to connect any of it together. I am a software developer, not a mathematician. So if anyone can walk me through some steps on how to get from the equation to 18 || 19, that would be great! 
This is a toy example representing a new Elliptic Curve Crypto operation where the actual modulus is $2^{256}$ large. So, trying all possible values x is not practical. WolframAlpha is capable of producing solutions to my large modulo equations in a fraction of a second so I know they aren't trying all possible values x.
Fermat’s Little Theorem seems the most promising so far, but I don't understand how to apply it to this equation. This post describes a solution but unfortunately their example is very basic and not very relatable to my equation.
Anything would be helpful here. Steps would be great. Thanks!

Let 
$$f(x)=x^4 -2x^3 + 5x^2 + 7x + 1\tag{1}$$
be defined over the finite field $\mathbb{F}_{23}$. Now check for a linear factor by checking for roots  over  $\mathbb{F}_{23}=\{0,\pm1,\pm2,\pm3,\pm4,\pm5\pm6,\pm7,\pm8,\pm9,\pm10,\pm11\}$. We find $f(-4)=f(-5)=0$, so $(x+4)$ and $(x+5)$ are linear factors. Now factor $f$ as two quadratics modulo $23$:
\begin{align*}
f(x)&=(x^2+9x-3)(x^2+ax+b)\\
&=x^4+(9+a)x^3+(9a-3+b)x^2+(9b-3a)x-3b
\end{align*}
Comparing the coefficients in $(1)$ for the powers of $x$:
\begin{array}\\
[x^3:] & -2=9+a\\
 [x^2:] & 5=9a-3+b\\
 [x:] & 7=9b-3a\\
 [const:]& 1=-3b\\
\end{array}
with $a$, $b$, $c$, $d\in\mathbb{F}_{23}$. Note this is a finite field so $-3b=1$ means $-3$ and $b$ are inverse mod $23$, making $b=15$. Now $a=-2-9=-11=12$ giving the factorization 
$$f(x)=(x^2+12x+15)(x+5)(x+4)$$
with the quadratic factor irreducible over $\mathbb{F}_{23}$ as it has no roots, since the discriminant of $(x^2+12x+15)$ is $15$ which is not a square modulo $23$.
