[F77] Metody numeryczne i problem z NaN

Witam. :)

Mam problem z poniższym kodem, który generuje (czym wzięta większa ilość elementów tablic (M) - tym "szybciej" się sypie) najpierw to co chcę, później cholendarnie wielkie liczby, na NaNach kończąc.

	DO i=1,M !M<=N

	!evolution
	DO tl=0,K
	call evol(psi_l(i),Hpsi_l(i))

	call evol(psi_l(i-1),Hpsi_l(i-1))
	call evol(psi_l(i+1),Hpsi_l(i+1))
	Hpsi_l(i)=Hpsi(psi_l(i),psi_l(i-1),psi_l(i+1),V_l)

	END DO !evolution over

	WRITE(*,*) psi_l(i)

	END

psi_l to uprzednio wygenerowany, stablicowany (tablica: psi_l(N), N=500) rozkład Gaussa N(0,4)
V_l to tablica stałych
Hpsi_l to stablicowana (wymiar tablicy oczywiście również =N);
Hpsi_l(i)=Hpsi(psi_l(i),psi_l(i-1),psi_l(i+1),V_l(i))
funkcja Hpsi:

	REAL  FUNCTION Hpsi(psi,psi1,psi2,V)
	REAL  psi,psi1,psi2,V

	Hpsi=-(psi1-2.*psi+psi2)+V*psi !druga pochodna psi ze wzoru trójpunktowego dla jednostkowej !odległości między węzłami +stała*psi

	END DO

Ewolucja:

	SUBROUTINE evol(psi,Hpsi)
	REAL psi,Hpsi

	psi=psi-Hpsi 

	END

Sypie się już dla K>=2.

Nie mam pojęcia dlaczego tak się dzieje - psi (Gauss) to małe wartości, Hpsi również (i przyjmuje poprawne wartości - tu nie ma źródła NaN), więc nie wiem skąd ta pętla z ewolucją (procedura ta dla pojedynczych elementów tablicy psi_l działa poprawnie) wyrzuca absurdalne liczby i NaN-y.

Gdzie może tkwić haczyk? Może należy jakoś inaczej sformułować tę główną pętlę?

Z góry dzięki za wszelkie porady. :)