for (int j = 0;j < N; j++){
t = j * dt;
k11 = f1(t,vx,x,vy,y);
k12 = f2(t,vx,x,vy,y);
k13 = f3(t,vx,x,vy,y);
k14 = f4(t,vx,x,vy,y);
k21 = f1(t+dt/2.0, vx+dt/2.0*k11, x+dt/2.0*k12, vy+dt/2.0*k13, y+dt/2.0*k14);
k22 = f2(t+dt/2.0, vx+dt/2.0*k11, x+dt/2.0*k12, vy+dt/2.0*k13, y+dt/2.0*k14);
k23 = f3(t+dt/2.0, vx+dt/2.0*k11, x+dt/2.0*k12, vy+dt/2.0*k13, y+dt/2.0*k14);
k24 = f4(t+dt/2.0, vx+dt/2.0*k11, x+dt/2.0*k12, vy+dt/2.0*k13, y+dt/2.0*k14);
k31 = f1(t+dt/2.0, vx+dt/2.0*k21, x+dt/2.0*k22, vy+dt/2.0*k23, y+dt/2.0*k24);
k32 = f2(t+dt/2.0, vx+dt/2.0*k21, x+dt/2.0*k22, vy+dt/2.0*k23, y+dt/2.0*k24);
k33 = f3(t+dt/2.0, vx+dt/2.0*k21, x+dt/2.0*k22, vy+dt/2.0*k23, y+dt/2.0*k24);
k34 = f4(t+dt/2.0, vx+dt/2.0*k21, x+dt/2.0*k22, vy+dt/2.0*k23, y+dt/2.0*k24);
k41 = f1(t+dt, vx+dt*k31, x+dt*k32, vy+dt*k33, y+dt*k34);
k42 = f2(t+dt, vx+dt*k31, x+dt*k32, vy+dt*k33, y+dt*k34);
k43 = f3(t+dt, vx+dt*k31, x+dt*k32, vy+dt*k33, y+dt*k34);
k44 = f4(t+dt, vx+dt*k31, x+dt*k32, vy+dt*k33, y+dt*k34);
vx = vx+dt/6.0*(k11+2.0*k21+2.0*k31+k41);
x = x+dt/6.0*(k12+2.0*k22+2.0*k32+k42);
vy = vy+dt/6.0*(k13+2.0*k23+2.0*k33+k43);
y = y+dt/6.0*(k14+2.0*k24+2.0*k34+k44);
repaint();