double f1 (double t, double x, double y,double step,double gamma,double omega ) { return step * ( -2 * gamma * y - pow(omega,2) * x) ; } double f2 (double t, double x, double y,double step,double gamma,double omega ) { return step * y ; } int main() { const double time =100; const double step =0.01; const unsigned int count = int((time)/step) + 1; vector time_vec(count); //vector new_x(count); //vector new_y(count); double X[(int)count]; double X1[(int)count]; double X2[(int)count]; double X3[(int)count]; double X4[(int)count]; double Y1[(int)count]; double Y2[(int)count]; double Y3[(int)count]; double Y4[(int)count]; double Y[(int)count]; double t[(int)count]; X[0]=5;Y[0] =5;t[0]=0; for(int i= 1 ; i < count ; i++) { double ntime = i * step ; X1[i]=f2(t[i-1],X[i-1],Y[i-1], step, gamma, omega); Y1[i]=f1(t[i-1],X[i-1],Y[i-1], step, gamma, omega); X2[i]=f2(t[i-1] + step/2, X[i-1]+ X1[i]/2.0,Y[i-1]+Y1[i]/2.0 ,step, gamma, omega); Y2[i]=f1(t[i-1] + step/2, X[i-1]+ X1[i]/2.0,Y[i-1]+Y1[i]/2.0 ,step, gamma, omega); X3[i]=f2(t[i-1] + step/2, X[i-1]+ X2[i]/2,Y[i-1]+Y2[i]/2,step, gamma, omega); Y3[i]=f1(t[i-1] + step/2, X[i-1]+X2[i]/2,Y[i-1]+Y2[i]/2,step, gamma, omega); X4[i]=f2(t[i-1] + step/2, X[i-1]+X3[i],Y[i-1]+Y3[i],step, gamma, omega); Y4[i]=f1(t[i-1] + step/2, X[i-1]+X3[i],Y[i-1]+Y3[i],step, gamma, omega); X[i]=X[i-1]+(X1[i]+2*X2[i]+2*X3[i]+X4[i])/6; Y[i]=Y[i-1]+(Y1[i]+2*Y2[i]+2*Y3[i]+Y4[i])/6; time_vec[i] = ntime; }