Неявный метод Рунге Кутта второго порядка (a)

// Ne9vn_RK_k2(a).cpp: îïðåäåëÿåò òî÷êó âõîäà äëÿ êîíñîëüíîãî ïðèëîæåíèÿ.
//

#include "stdafx.h"
#include <stdlib.h>
#include <math.h>

struct y1y2
{
    double y1, y2, h;
};

double f1(y1y2 y1_y2)
{
    return (-501*y1_y2.y1+500*y1_y2.y2);
}

double f2(y1y2 y1_y2)
{
    return (500*y1_y2.y1-501*y1_y2.y2);
}

void print(int i, y1y2 y1_y2)
{
    printf("hod %d y1=%f y2=%f\n",i,y1_y2.y1,y1_y2.y2);
    system("pause");
}

y1y2 eler(y1y2 y1_y2_h)
{
    y1y2 yn;
    yn.y1=y1_y2_h.y1+y1_y2_h.h*f1(y1_y2_h);
    yn.y2=y1_y2_h.y2+y1_y2_h.h*f2(y1_y2_h);
    yn.h=y1_y2_h.h;
    return yn;
}

y1y2 pr_it(y1y2 y1_y2_h, y1y2 y1ny2nh)
{
    y1y2 n={0, 0, 0}, t={0, 0, 0}, y1ky2kh=y1ny2nh;
    double eps=0.0001;
    do
    {
        y1ny2nh=y1ky2kh;
        n.y1=f1(y1ny2nh);
        n.y2=f2(y1ny2nh);
        t.y1=y1ny2nh.y1-y1_y2_h.h*n.y1;
        t.y2=y1ny2nh.y2-y1_y2_h.h*n.y2;
        y1ky2kh.y1=y1_y2_h.y1+(y1_y2_h.h/2)*(n.y1+f1(t));
        y1ky2kh.y2=y1_y2_h.y2+(y1_y2_h.h/2)*(n.y2+f2(t));
    }
    while((fabs(y1ky2kh.y1-y1ny2nh.y1)>eps)||(fabs(y1ky2kh.y2-y1ny2nh.y2)>eps));
    return y1ky2kh;
}

y1y2 RK(y1y2 y1_y2_h, y1y2 y1ky2kh)
{
    y1y2 n={f1(y1ky2kh),f2(y1ky2kh),0}, t={y1ky2kh.y1-y1_y2_h.h*n.y1,y1ky2kh.y2-y1_y2_h.h*n.y2,0};
    y1ky2kh.y1=y1_y2_h.y1+(y1_y2_h.h/2)*(n.y1+f1(t));
    y1ky2kh.y2=y1_y2_h.y2+(y1_y2_h.h/2)*(n.y2+f2(t));
    return y1ky2kh;
}

void main(void)
{
    y1y2 y1_y2_h={5, -10, 0.000001};
    int i=0;
    print(i, y1_y2_h);
    while(1)
    {
        y1_y2_h=RK(y1_y2_h, pr_it(y1_y2_h, eler(y1_y2_h)));
        print(((++i)*10), y1_y2_h);
    }
}