Написать программу для вычисления волновых функций Кулона на Fortran 90

Необходимо написать программу, которая бы вычисляла волновые функции водородоподобного атома в случае непрерывного спектра. У меня есть определенные наработки. Но программа не всегда дает верный результат. Можно попытаться улучшить предложенный мною алгоритм, либо использовать другой.

В pdf файлах "shablon" и "thpu" описан алгоритм и приведена функция на Fortran 90 для вычисления. Но она работает не всегда верно:

В обозначениях файла "shablon". При l=1 программа работает правильно, хотя иногда не тот знак получается в результате. А при l>1 программа выдает верный результат только при больших значениях радиального расстояния r.

Также прикреплю проект "my_project1", в котором будут проведены некоторые тестовые вычисления. В файле проекта "coulomb_thpu.f90" содержится функция, описанная в pdf "thpu". Файл "coulomb_barnet.for" не нужен, но может быть использован для проверки.

В файле "peng2010". Будут таблицы с вычисленными волновыми функциями при l=1,10,80,200. Эти данные можно использовать для проверки вычислений.