Да и ряд

не шибко-то сходится.
Ну это, судя по всему, ошибка. Ряд и не должен ВЕЗДЕ сходиться. В точке, где источник, он просто обязан расходиться. Там же должна получиться бесконечность. Ну и пусть ТАМ расходится, это правильно. В остальных точках осциллирующие множители приведут к тому, что сходимость все же будет.
Кстати, хорошо бы избавиться от этой плохой сходимости вблизи источника (и, как следствие, необходимости брать очень много членов ряда в практических вычислениях). Но с этим в третьем методе будем бороться (а есть еще и четвертый, хоть и не метод решения, но метод улучшения ряда). А пока второй. Будем избавляться от ДВОЙНОГО ряда, делать одинарный.
Собственно, нужно просто вернуться к первоначальным рассуждениям ТС, только немного их модифицировать. Там возникла проблема с функцией Неймана. Почему она возникла? Дело в том, что обычное в цилиндрических задачах условие регулярности в нуле --- это тоже, своего рода, граничное условие. А тут сингулярность источника "наехала" на это граничное условие. Особый случай! В принципе, можно попытаться "пробить", рассматривая асимтотику. Но лучше бы от этого всего избавиться. Как? А очень просто: разложить не по

, а потом решать уравнение Бесселя (неоднородное), а сделать все наоборот: разложить по круглым (обычным) бесселям

и решать неоднородное уравнение вида:

Здесь сингулярность источника уже не "наезжает" на граничное условие.
Я написал для заряда, а не для диполя (и опустил, индексы и множители). Но все абсолютно аналогично, переделать на диполь --- самостоятельное упражнение.
Как решать уравнение (*)? Можно, конечно, методом вариации постоянных (или использовать, в качестве частного решения, решение на всей оси, получаемое методом интеграла Фурье), но лучше (и проще) иначе. Сейчас опишу простой прием, позволяющий быстро получить решение, причем в очень "хорошей" форме.
Заметим, что везде КРОМЕ ТОЧКИ

уравнение является однородным, два базисных решения

которого известны (находятся на раз). Поэтому можно написать:

Где

--- функция Хевисайда (единичный скачок). Подставляем этот анзац в (*) и вспоминаем, что

. Тогда в левой части получится сумма двух сингулярных членов (с коэффициентами):

и

. Чтобы уравнение выполнялось, нужно чтобы выполнялось два (линейных!) условия на коэффициенты.
Два условия на четыре коэффициента

есть. Еще два граничных условия, и все

определяются запросто!
В качестве дополнительного упражнения (и проверки) рекомендую стандартным образом разложить полученную функцию

в ряд Фурье, и убедиться, что получается в точности то, что ранее было получено методом разложения по собственным функциям лапласиана.