Метод Рунге-Кутты системы ДУ

Нужно написать программу для решения задачи Коши для системы дифференциальных уравнений типа dy/dt=A11*y+A12*x+B1 dx/dt=A21*y+B22*x+B2 И почему то получается линейная зависимость. Не правильно написана программа или может алгоритм?

#include <iostream>
#include "conio.h"
#include <stdio.h>
#include <math.h>
#include <fstream>
using namespace std;
double f2(double x, double y,double t)
{   
    return  (-53.3333*y-7333.3333*x+40000);
}
double f1(double x, double y,double t)
{
    return  (73.333*y-3666.6667*x+2000);
}
int main()
{   
ofstream outfile("out.txt");
ofstream outfile1("out1.txt");
ofstream outfile2("out2.txt");
   system("mode con cols=126 lines=1020");
    int i;
    double x, y, d, k, k1y,k2y,k3y,k4y, dy,dx,dz,k1x,k2x,k3x,k4x,t=0;
    double x0=-0;x=x0;
   double y0=0;y=y0;
   double h=0.0239/1000;
    for(i=0;i++,i<=1000;t=t+h,y=y+dy,x=x+dx)
    {       
            k1x=h*f1(x,y,t);              
            k1y=h*f2(x,y,t);
            k2x=h*f1(x+k1x/2,y+k1y/2,t+h/2);
            k2y=h*f2(x+k1x/2,y+k1y/2,t+h/2);
            k3x=h*f1(x+k2x/2,y+k2y/2,t+h/2);
            k3y=h*f2(x+k2x/2,y+k2y/2,t+h/2);
            k4x=h*f1(x+k3x,y+k3y,t+h/2);
            k4y=h*f2(x+k3x,y+k3y,t+h/2);
            dy=(k1y+2*k2y+2*k3y+k4y)/6;
            dx=(k1x+2*k2x+2*k3x+k4x)/6;
            cout<<"x="<<x<<"   "<<"y="<<y<<"    ""t="<<t<<"\n";
            outfile <<x<<";"<<y<<";"<<t<< endl;
            //outfile <<"{"<<x<<","<<y<<","<<"}"<< endl;
             outfile1 <<y<<";"<< endl;
             outfile2 <<x<<";"<< endl;
    }
   return 0;
}