Как решить систему ОДУ методом Рунге-Кутты 4 порядка? [требует правки]

Привет, хочу решить систему ОДУ методом Рунге-Кутты 4 порядка. Правильно ли реализован этот метод?

public class RealizationMethods implements Methods{
double a = 0, b = 5;
int n = 10, k;
double k11, k21, k31, k41, k12, k22, k32, k42, k13, k23, k33, k43, k14, k24, k34, k44, k15, k25, k35, k45;
double h, y1, y2, y3, y4, y5, yn1, yn2, yn3, yn4, yn5, fi = 0, dfi = 0, dzeta = 1;
@Override
public double f1(double etta, double y2) {
    return y2;
}
@Override
public double f2(double etta, double y3) {
    return y3;
}
@Override
public double f3(double etta, double y1, double y2, double y3) {
    return -y1*y3+y2*y2-1;
}
@Override
public double f4(double etta, double y5) {
    return y5;
}
@Override
public double f5(double etta, double y1, double y2, double y4, double y5) {
    return -(y1*y5-y2*y4);
}
@Override
public double RK() {
    double etta;
    h = (b - a)/n;
    for(int i = 0; i < n-1; i++) {
        etta = a+i*h;
        k11 = h * f1(etta, y2);
        k21 = h * f1(etta + h / 2.0, y2 + k11 * h / 2.0);
        k31 = h * f1(etta + h / 2.0, y2 + k21 * h / 2.0);
        k41 = h * f1(etta + h, y2 + h*k31);
        k12 = h * f2(etta, y3);
        k22 = h * f2(etta + h / 2.0, y3 + k12 * h / 2.0);
        k32 = h * f2(etta + h / 2.0, y3 + k22 * h / 2.0);
        k42 = h * f2(etta + h, y3 + h * k32);
        k13 = h * f3(etta, y1, y2, y3);
        k23 = h * f3(etta + h / 2.0, y1 + k13 * h / 2.0, y2 + k13 * h / 2.0, y3 + k13 * h / 2.0);
        k33 = h * f3(etta + h / 2.0, y1 + k23 * h / 2.0, y2 + k23 * h / 2.0, y3 + k23 * h / 2.0);
        k43 = h * f3(etta + h, y1 + h * k33, y2 + h * k33, y3 + h * k33);
        k14 = h * f4(etta, y5);
        k24 = h * f4(etta + h / 2.0, y5 + k14 * h / 2.0);
        k34 = h * f4(etta + h / 2.0, y5 + k24 * h / 2.0);
        k44 = h * f4(etta + h, y5 + h * k34);
        k15 = h * f5(etta, y1, y2, y4, y5);
        k25 = h * f5(etta + h / 2.0, y1 + k15 * h / 2.0, y2 + k15 * h / 2.0, y4 + k15 * h / 2.0, y5 + k15 * h / 2.0);
        k35 = h * f5(etta + h / 2.0, y1 + k25 * h / 2.0, y2 + k25 * h / 2.0, y4 + k25 * h / 2.0, y5 + k25 * h / 2.0);
        k45 = h * f5(etta + h, y1 + h * k35, y2 + h * k35, y4 + h * k35, y5 + h * k35);
        yn1 = y1 + h/6*(k11 + 2*k21 + 2*k31 + k41);
        yn2 = y2 + h/6*(k12 + 2*k22 + 2*k32 + k42);
        yn3 = y3 + h/6*(k13 + 2*k23 + 2*k33 + k43);
        yn4 = y4 + h/6*(k14 + 2*k24 + 2*k34 + k44);
        yn5 = y5 + h/6*(k15 + 2*k25 + 2*k35 + k45);}
        System.out.println(yn1);
        System.out.println(yn2);
        System.out.println(yn3);
        System.out.println(yn4);
        System.out.println(yn5);
    fi = y1;
    dfi = y2;
    dzeta = y4;
    return 0;
}