I have a system of 4 nonlinear first-order ODEs: $$ \begin{align} \frac{dS}{dt} &= -\beta SI \\ \frac{dE}{dt} &= \beta SI - \sigma E \\ \frac{dI}{dt} &= \sigma E - \gamma I \\ \frac{dR}{dt} &= \gamma I, \end{align} $$ where $\beta$, $\sigma$ and $\gamma$ are positive parameters of the system.
I found that the steady-state (fixed point) will be a line that is defined by $I=0$, $E=0$ (considering only 3D $S-E-I$ space since $N=S+E+I+R$ remains constant).
I constructed the Jacobian matrix: $$ \begin{align} F(S, I)&= -\beta SI \\ G(S, E, I)&= \beta SI - \sigma E \\ H(I, E)&= \sigma E - \gamma I \\ J(F, G, H)= \begin{pmatrix} \frac{\partial F}{\partial S} & \frac{\partial F}{\partial E} & \frac{\partial F}{\partial I} \\ \frac{\partial G}{\partial S} & \frac{\partial G}{\partial E} & \frac{\partial G}{\partial I} \\ \frac{\partial H}{\partial S} & \frac{\partial H}{\partial E} & \frac{\partial H}{\partial I} \end{pmatrix} &= \begin{pmatrix} -\beta I & 0 & -\beta S \\ \beta I & -\sigma & \beta S \\ 0 & \sigma & -\gamma \end{pmatrix} \end{align} $$
This matrix is evalueted at the fixed line: $$ J(S, 0, 0)= \begin{pmatrix} 0 & 0 & -\beta S \\ 0 & -\sigma & \beta S \\ 0 & \sigma & -\gamma \end{pmatrix} $$
This returns 3 eigenvalues: $$ \begin{align} \lambda_1&= 0 \\ \lambda_{2,3}&=\frac{\sigma +\gamma\pm\sqrt{\sigma^2+2\sigma\gamma+ \gamma^2+4\sigma(\beta S-\gamma)}}{2} \end{align} $$
How do I interpret these eigenvalues in order to determine the type of point (sink, source, saddle, etc)?
And bonus question: Is there any way to make a nice 3D phaseplot? For 2D systems, I used matplotlib's streamplot function.
Since zero is an eigenvalue (as it will always be when you have a whole line of fixed points), the standard linearization theorem (Hartman–Grobman) unfortunately doesn't tell you anything about the type. For that theorem to be applicable, the equilibrium needs to be hyperbolic, meaning that there are no eigenvalues with real part zero. But see the comment by Ian below for more advanced methods that may sometimes be helpful.
However, regarding this particular problem, a better way to understand the system is to use the constant of motion $S+E+I-\frac{\gamma}{\beta}\,\ln S$, which basically plays the same role as the constant of motion $S+I-\frac{\gamma}{\beta}\,\ln S$ does in the SIR model ($dS/dt=-\beta SI$, $dI/st=\beta SI-\gamma I$, $dR/dt=\gamma I$), if you're familiar with that (which you probably should be before you study the SEIR model).