Перейдите к вариационной постановке задачи:

Для граничных условий третьего типа

выпишите квадратичную форму оператора



Теперь сажайте всё это на сетку, минимизируйте полученную квадратичную форму -- и получите замкнутую систему разностных уравнений. (Граничные условия второго и третьего типа, в отличие от первого, являются
следствием минимизации квадратичной формы, на область определения которой изначально никаких граничных условий не накладывалось.)
(если ничего не напутал в знаках; ну разберётесь, если что)