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

[DenProx]

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

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

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

Очень нужна ваша помощь!

Есть уже некоторые наработки:

void runge_k(
    void f(double *y, double *ys, double t),
    double *y, int n,
    duble tn, double tk,
    int m, double delt)
{
    double k1[n],k2[n],k3[n],k4[n],Yh[n+1];
    double h=delt;
    double* ys=y+n+1;
    double X;
    for(X=tn;X<tk;X+=h){
        f(y,ys,X);
        
        for( int i=0;i<n;i++) k1[i]=h*y[i+1];
        for( int i=0;i<n;i++) Yh[i]=y[i]+0.5*k1[i];         
        f(Yh,Yh+n+1,X+h*0.5);
        
        for( int i=0;i<n;i++) k2[i]=h*Yh[i+1];
        for( int i=0;i<n;i++) Yh[i]=y[i]+0.5*k2[i];
        f(Yh,Yh+n+1,X+h*0.5);
        
        for( int i=0;i<n;i++) k3[i]=h*Yh[i+1];
        for( int i=0;i<n;i++) Yh[i]=y[i]+k3[i];
        f(Yh,Yh+n+1,X+h);
        
        for( int i=0;i<n;i++) k4[i]=h*Yh[i+1];
        
        for( int i=0;i<n;i++) y[i]+=1./6.*(k1[i]+2*k2[i]+2*k3[i]+k4[i]);
    }
    h=tk-X;
    f(y,ys,X);
    
    for( int i=0;i<n;i++) k1[i]=h*y[i+1];
    for( int i=0;i<n;i++) Yh[i]=y[i]+0.5*k1[i];         
    f(Yh,Yh+n+1,X+h*0.5);
    
    for( int i=0;i<n;i++) k2[i]=h*Yh[i+1];
    for( int i=0;i<n;i++) Yh[i]=y[i]+0.5*k2[i];
    f(Yh,Yh+n+1,X+h*0.5);
    
    for( int i=0;i<n;i++) k3[i]=h*Yh[i+1];
    for( int i=0;i<n;i++) Yh[i]=y[i]+k3[i];
    f(Yh,Yh+n+1,X+h);
    
    for( int i=0;i<n;i++) k4[i]=h*Yh[i+1];
    
    for( int i=0;i<n;i++) y[i]+=1./6.*(k1[i]+2*k2[i]+2*k3[i]+k4[i]);
}

[keiya]

программу проверили? нужно взять какое нибудь известное точное решение и сравнить с численным