Применить теорему Гаусса ко всей сферической оболочке?
Да, верно.
В чем же все-таки дело?
Пара предварительных замечаний.
1. Как Вы знаете, потенциал

определён с точностью до константы. Выбор того или иного потенциала из бесконечного множества функций, отличающихся константой, называется калибровкой. Пусть система зарядов финитна (вне некоторого шара зарядов нет), но суммарный заряд

ненулевой. Тогда, очевидно, изменение калибровки

приведёт к изменению значения интеграла:

ЛЛ при выводе существенно использовали то, что на бесконечности потенциал стремится к нулю — там, где отбрасывали объёмный интеграл с дивергенцией, сводящийся к интегралу по бесконечно удалённой поверхности. Так что формулу

можно использовать лишь при

, что и определяет правильную калибровку. В нашей задаче это не так.
Однако

парадокс этим не снимается, так как в нашем случае

, и калибровка не должна влиять на результат.
2. Равенство

гарантируется только для интегралов по всему пространству, но для произвольной области интегрирования значения левой и правой части будут различны. Можно сказать, что «с точки зрения» первого и второго выражения электростатическая энергия локализована в пространстве по-разному. Вопрос об истинной локализации энергии (и её существовании) сложный и не должен сейчас обсуждаться. Сейчас предлагается просто считать, что оба интеграла имеют физический смысл

, лишь когда применяются ко всему пространству. Применять их к части пространства можно только в математическом смысле. Фраза
В этой области столько-то энергии с точки зрения выражения 
будет относиться лишь к аспекту вычисления.
Так вот, с точки зрения выражения

существенная часть энергии содержится на поверхностях. Об этом легко забыть. Коль вспомнили, как её учесть? Можно применить аппарат обобщённых функций, упомянутый
Munin, тогда в выражение для

на поверхностях будет входить дельта-функция. А можно догадаться, к чему сведётся дополнительное слагаемое «на человеческом языке»:

В нашем случае на внутренней поверхности потенциал нулевой. Зато «энергия» на внешней поверхности по модулю превосходит «энергию» непрерывного распределения и в сумме с ней даёт правильный результат. Вычисляется интеграл легко, так как на поверхности и

, и

константы.