I will follow the original paper as much as I can. So, let us consider the solution of the equation

$$u_t-6uu_x+u_{xxx}=\epsilon u.$$

You will have a general solution $u=u(x,t;\epsilon)=\sum_{n=0}^\infty\epsilon^n u_n(x,t)$ after a power expansion on $\epsilon$. Now, turning your attention to the Schroedinger equation

$$\psi_{xx}-[u(x,t,;\epsilon)-\lambda]\psi=0$$

you will notice that the problem is now

$$\psi_{xx}-[u_0(x,t)+\epsilon u_1(x,t)+\ldots-\lambda]\psi=0$$

and this is amenable to standard Rayleigh-Schroedinger perturbation theory. This means that the first correction preserves the discrete nature of the spectrum and will be real yet. You can iterate this procedure going to higher orders in $\epsilon$ still preserving the inverse scattering method. This can be seen in the follwing way. From the given Schroedinger equation one has

$$u(x,t;\epsilon)=\lambda+\frac{\psi_{xx}}{\psi}.$$

The condition granting that this solves KdV equation is $\lambda_t=0$ when $\psi\rightarrow 0$ given $|x|\rightarrow\infty$ and $u$ evolves by KdV equation. This must be true also for the approximate solution. So, from this equation we can evaluate the next-to-leading order correction by setting

$$\lambda=\lambda_0+\epsilon\lambda_1+\ldots$$
$$\psi = \psi_0+\epsilon\psi_1+\ldots$$

and we get

$$u(x,t;\epsilon)=\lambda_0+\epsilon\lambda_1+\frac{\psi_{0xx}}{\psi_0}-\epsilon\frac{\psi_1\psi_{0xx}-\psi_0\psi_{1xx}}{\psi_0^2}+\ldots$$

and so

$$u_0(x,t)=\lambda_0+\frac{\psi_{0xx}}{\psi_0}$$

$$u_1(x,t)=\lambda_1-\frac{\psi_1\psi_{0xx}-\psi_0\psi_{1xx}}{\psi_0^2}$$

$$\vdots$$

Turning back to KdV equation for the leading order one has

$$u_{0t}-6u_0u_{0x}+u_{0xxx}=0$$

and this is exactly the problem solved by inverse scattering method that we recover as it should. The next-to-leading order correction gives:

$$u_{1t}+6(u_0u_{1x}+u_1u_{0x})+u_{1xxx}=u_0$$

that is a linear equation. We notice here that we are consistent with the condition $\psi_1\rightarrow 0$ when $|x|\rightarrow\infty$ provided this is true for $\psi_0$ and this grants also $u_1\rightarrow 0$ as it should. One can check at this point that this is equivalent on doing the same perturbation approach on the scattering problem for the inverse scattering method (see eq.(5) in the original paper by Miura et al.).