Итак, приступим. Рассматриваем уравнение вида

Пусть

--- наименьшее решение ассоциированного уравнения Пелля

.
Теорема. Базовые решения

уравнения

определяются условиями:

в случае

; если же

, то

а также

Базовость означает, что

, где

(обычно берут

, но при указанном значении

интервал для

получается самым узким). Произвольное решение уравнения

представляется в виде

, где

, а

--- одно из базовых. В качестве примера решим уравнение

где

(

произвольно,

фиксировано). Имеем

, поэтому для

имеем оценку

, где

не зависит от

и

(но зависит от

). Фиксируем одно из таких

. Тогда для

имеем

, где

зависит только от

. Значит,

, где

. Подставляем это в уравнение и получаем

Таким образом, имеется конечное число линейных (относительно

) выражений для

. Оставим из них те, для которых

при всех достаточно больших

(нижняя граница для

зависит только от

; маленькие

разбираем отдельно). Для каждого из оставшихся выражений

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

(в роли

--- константы

, в роли

--- линейные выражения относительно

) и пишем ответ. (Небрежность в оценках здесь есть, так что проверяйте.)
Теперь наша задача (вот тут я старался быть аккуратным). Она сводится к исследованию уравнения

где

,

,

, при этом

. Здесь

При ограничениях на

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

двойное неравенство

Но тогда для

с учётом ограничений на

мы получили бы

откуда

и, как следствие, уравнение

Дискриминант этого уравнения должен быть точным квадратом:

Вновь учитывая ограничения на

, приходим к неравенствам

Значит,

и

. Таким образом,

Все искомые пары

натуральных чисел теперь можно найти из формулы

где

(и, конечно,

--- произвольно).
Не вижу, почему бы всё это не прошло и для других интервалов для

.