imaginary time propagation to find ground state wavefunction

I understand the basic idea of imaginary time propagation method:

The wavefunction $\psi(x,t)$ as a superposition of energy eigenstates $\phi_m(x)$: $$\psi(x,t)=\sum_m \phi_m(x)e^{-iE_mt/\hbar}$$ In imaginary time, ${t}\rightarrow{-it}$, $$\psi(x,-it)=\sum_m \phi_m(x)e^{-E_mt/\hbar}=e^{-E_ot/\hbar}\bigg(\phi_o+\phi_1e^{-(E_1-E_o)t/\hbar}+\phi_1e^{-(E_2-E_o)t/\hbar}+........\bigg)$$ As time $t$ increases the terms with greater exponents will decay more rapidly, leaving behind the ground state.

But how do i make a program or subroutine to find the ground state wavefunction using imaginary time propagation for a basic one dimensional problem, preferably in Fortran ?

Ref to 7.12 in an introduction to computational physics by Pang.

The macroscopic state of Bose–Einstein condensate is described by the Gross–Pitaevskii equation: $$i\hbar\frac{\partial \Psi(x,t)}{\partial t}=\frac{-\hbar^{2}}{2m}\frac{\partial^2 \Psi(x,t)}{\partial x^2}+V_{ext}\Psi(x,t)+g|\Psi(x,t)|^2\Psi(x,t)$$ where $\psi(x,t)$ is the macroscopic wavefunction, $m$ and $g$ are two system-related parameters, and $V_{ext}(x)$ is the external potential. Develop a program that solves this equation numerically. Consider various potentials, such as a parabolic or a square well.

Can I apply Finite-difference method or Crank-Nickelson method, similar to simple diffusion equation, by substituting for $\frac{\partial \Psi(x,t)}{\partial t}, \frac{\partial^2 \Psi(x,t)}{\partial x^2}, |\Psi(x,t)|^2\Psi(x,t), \Psi(x,t)$, after inputting the imaginary time to get rid of '$i$' ?

or how do I apply imaginary time evolution to Schrodinger equation of very basic 1D system ?

Note: It'd be very helpful if someone can give me advice on where to start.

I believe -- it's been a while since I worked on it -- that you have to move to imaginary time (as you said, by using $-it$ instead of $t$). Then you want to propagate along time the equation obtained (something like : $\frac{\partial \Psi}{\partial t} = \frac{1}{2} \Delta \Psi - V(x,y)\Psi - \beta |\Psi|^2 \Psi$ ), which can be done as follows.
You can then solve these "sub-timesteps" on $|\Psi|^2 = \rho$. One is easier to solve in Fourier space, the other is exactly solvable.