Нужно составить программу для решения дифф. уравнений методом Рунге-Кутта

Сами как думаете - кто Вам тут что напишет? Я, конечно же, только самое простое:

 from math import* 

fun, a, b, ya, ε = input('f(x,y) a b y(a) ε: ').split() 

f, a, b = lambda x, y: eval(fun), float(a), float(b) 

n, Y = int((b - a) / float(sqrt(float(ε)))), [float(ya)] 

h = (b - a) / n; hp = h * 0.5 

for i in range(n): 

    y = Y[i]; x = a + i * h 

    Y.append(y+h*f(x+hp, y+hp*f(x,y))) 

for i in range(n + 1): 

	x = a + i * h 

	print('%7.4f%22.15f' % (x, Y[i])) 

Это с методом Рунге-Кутты второго порядка. Для больших интервалов надо, естественно, что-нибудь поточнее. Вышеприведённый код легко модифицировать, доведя метод РК до хотя бы четвёртого порядка, но я думаю, что это уже излишество! Тут решение получается просто в виде таблицы узловых значений функции. Все необходимые исходные данные вводятся с экрана через пробел (а сама функция f(x,y) должна быть без пробелов !), поэтому можно проверить точность, например, для экспоненты на интервале [0;1], введя:
y 0 1 1 1e-6
(T.e. решить ОДУ у'=у, у(0)=1 на [0;1] с точностью 1е-6)
Абсолютная ошибка будет меньше 1е-6.
А вообще надо указывать что конкретно нужно: например, порядковую точность метода, а не ε-точность и, может быть, не требование получения списка узловых значений, а задание на получение работающей функции с узловыми значениями в узлах и восполненную сплайнами по её узловым значениям в междуузлиях. И ещё, может быть, надо визуализировать полученные данные. Но графика - это вообще особая тема! Кстати, а как там насчёт библиотек? Что - разве в Пайтоне нет уже готовых инструментов для решения такой вот ерунды? Да их там полным-полно!