Неявный метод Рунге Кутта второго порядка (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);
        }
}