Численное интегрирование системы дифференциальных уравнений методом Рунге - Кутта - вопрос №408462

Доброго времени суток. У меня есть вот такая задачка:

Разработать функцию для численного интегрирования системы дифференциальных уравнений методом Рунге — Кутта. Прототип функции:

void runge_k(void f(double *y, double *ys, double t), double *y, int n, duble tn, double tk, int m, double delt);

где:

f — функция вычисления правых частей системы дифференциальных уравнений;
y — массив размера n значений зависимых переменных;
ys — массив размера nзначений производных;
n — порядок системы дифференциальных уравнений;
t — независимая переменная;
tn — начальное значение интервала интегрирования;
tk — конечное значение интервала интегрирования;
m — начальное число разбиений отрезка интегрирования ;
delt — шаг интегрирования.

Шаг интегрирования для метода использовать 0,0001.

необходимо применить эту функцию для интегрирования дифференциального уравнения 3-го порядка:

y''' + 2y'' + 3y' + y = 5 + x^2

в интервале x [0,2] с шагом delta X = 0,1
и начальными условиями: x=0; y(0) =1; y'(0) = 0,1; y''(0) = 0

Смысл в том, что ни один программист которого я знаю, не взялся ее решить. И ни на одном форуме, на которых я пытался найти человека который смог бы решить данную задачу, я не нашел. 

Если кто то в состоянии ее решить, для интереса (не за деньги), буду только рад :) 

p.s. решение интересует на Си или Pascal