Linearizing the Nonlinear Schrodinger equation near Stationary Solutions Problem Statement:
Given the Lugiato-Lefever equation, linearize the equation and determine the dynamics near a stationary solution by looking for a stationary solution with a small perturbation.
Relevent Equations:
The  Lugiato-Lefever Equation is given by:
$$ \frac {d\psi} {dt} = -(1 + i \alpha)\psi + \frac {i} {2} \frac {d^2\psi} {dx^2} + i \psi^2 \bar \psi + F $$ 
It can be shown that when looking for stationary solutions in both time and space, the following is a solution:
$$ \psi_e = \frac F  {1 + i (\alpha - \rho)} $$
When $\rho$ is the solution of the polynomial:
$$ F ^ 2 = (1 + (\rho - \alpha) ^ 2)\rho $$
Attempt at a solution:
I started by substituting the following anzatz:
$$ \psi = \psi_e + \psi_1  $$
When $ |\psi_1| \ll 1 $. Substituting the above into the equation yields:
$$ \frac {d\psi_1} {dt} = -(1 + i\alpha)\psi_1 + \frac i 2 \frac {\partial ^ 2 \psi _1 } {\partial x ^ 2 }  + i (\bar \psi_1 \psi_1 ^2 + \bar \psi_e \psi _1 ^2 + 2\bar\psi_e \psi_e \psi_1 + \bar \psi_1 \psi_e ^2 + 2\bar\psi_1 \psi_e \psi_1) $$
I know $\psi_e$ as it is the stationary solution from above, and for $\psi_1$ I can look for solutions of the form:
$$ \psi_1 = e^{i\omega t}(a e^{ikx} + \bar b e^{-ikx})$$
I am not sure if I am in the right direction because I was told at this point I should only take the first order in the perturbation ($\psi_1$), yet I am supposed to get a coupled system of equations on a and b that depends on k, and this will only happen with the terms quadratic in $ \psi_1 $. Another thing that worries me is that usually in treatments like this I am only supposed to get an equation on $\psi_1 $ yet I still have terms dependent on $\psi_e $.

I think your equation is formulated incorrectly. The right Lugiato-Lefever equation takes the form (see here) 
$$
\frac{\partial\psi}{\partial t}=-i\theta\psi+\frac{i}{2}\frac{\partial^2\psi}{\partial x^2}+i|\psi|^2\psi+F.
$$
Please, note the use of the partial derivatives and the module in the non-linear term.
Finally, consider to use a signpost $\epsilon$ for your perturbed solution and write
$$
\psi=\psi_e+\epsilon\psi_1+\ldots.
$$
To get the right first order equation you just need to omit all the powers of $\epsilon^2$ and higher. At the end of computation just put $\epsilon=1$ and you are done.
Can you take from here?
