Как мат. пакеты вычисляют интегралы?

Гугли "Компьютерная алгебра".
Но алгоритмы там недоступны для понимания большинству даже программистов.

Не объяснишь профану - что такое "символьные методы решения" хоть чего-нибудь?

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

Очень примерный пример:

Пусть у нас есть ∫3+x^2*9dx

Для начала будет использовано свойство суммы:

∫3dx+∫x^2*9dx

Затем 2-ой интеграл приведен к стандартному виду:
∫3dx+∫9x^2dx

Еще одно свойство:

∫3dx+9∫x^2dx

Теперь будут выбраны шаблоны и произведены действия:

3x + const + 9*(x^3/3) + const

Затем будет произведено упрощение полученного выражения:

3x + 3x^3 + const

Затем его приведение к стандартному виду:

3x^3 + 3x + const

Затем, возможно, дополнительное упрощение:

3x(x^2+1) + const

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

S{x=0;+∞}x^(-a•x)dx, a>0

S{x⁴exp[√sin(x²)]}dx