Построены вычислительные схемы решения задачи Штурма-Лиувилля с однородными краевыми условиями первого, второго и третьего рода методом конечных элементов, сохраняющие в приближённых решениях свойства непрерывности производных искомых решений. Выведены рекуррентные соотношения для вычисления в аналитическом виде интерполяционных полиномов Эрмита с узлами произвольной кратности. Из интерполяционных полиномов Эрмита сконструированы базисные кусочно–полиномиальные функции на конечноэлементной сетке с переменным шагом, аппроксимирующие решение исходной задачи. Исходная задача Штурма-Лиувилля в базисе кусочно полиномиальных функций редуцируется к обобщённой алгебраической задаче на собственные значения с ленточными матрицами жёсткости и масс. Построены матрицы жёсткости и масс в виде сумм интегралов, содержащих заданные коэффициентные и потенциальные функции исходного самосопряжённого дифференциального уравнения и вычисленные интерполяционные полиномы Эрмита и их производные. Интегрирование выполняется с помощью гауссовых квадратур, а в специальных случаях, включающих кусочно-полиномиальные коэффициентные и потенциальные функции, в аналитическом виде. Эффективность и скорость сходимости предложенных вычислительных схем и разработанных алгоритмов и программ в среде Maple-Fortran доказана численным анализом тестовых расчётов точно решаемых задач Штурма-Лиувилля с непрерывными и кусочно-непрерывными потенциальными функциями.