И описывать сеточные уравнения для именно этих узлов отдельно, сшивая там производные (трёхточечные разностные).
обьясните пожалуйста вот сдесь поподробнее
Стандартно уравнение в каждом узле аппроксимируется просто заменой каждой производной на соотв. конечноразностную. А если

-- точка разрыва, то для этого узла надо не уравнение аппроксимировать, а приравнять друг дружке её левую и правую производные (по узлам, соответственно,

и

.
Правда, система при этом получится формально не трёх-, а пятидиагональная. Но, во-первых, выпирающие ненулевые элементы для каждого узла легко убить вручную соответствующим шагом метода Гаусса. А во-вторых, если лень думать, то можно легко обобщить процедуру, реализующую метод прогонки, на пятидиагональные матрицы. Правда, при этом время счёта в несколько раз увеличится, но не на порядок же; так что если время не критично, то какая разница.