>>108
起き抜けに作ってみた。とりあえずこれ参考にしてくれ。

#define P_X 0
#define P_Y 1
#define P_VX 2
#define P_VY 3
void calc(double kn[4],double p[4]){
double r = sqrt (p[P_X] * p[P_X] + p[P_Y] * p[P_Y] );
kn[P_X] = p[P_VX];
kn[P_Y] = p[P_VY];
kn[P_VX] = - p[P_X] / (r * r * r);
kn[P_VY] = - p[P_Y] / (r * r * r);
}
int main(int argc, char * argv[]){
double w[4],p[4],k1[4],k2[4],k3[4],k4[4];
int iter = 1000;
double T = 50;
double dt = T / iter;
int t,i; /* t 計算中の時刻 */
w[P_X] = 2.0; w[P_Y] = 0.0; w[P_VX] = 0.0; w[P_VY] = 0.5;
for( t = 0 ; t < iter ; t++){
for(i=0;i<4;i++){ p[i] = w[i] ; } calc(k1,p);
for(i=0;i<4;i++){ p[i] = w[i] + 0.5 * dt * k1[i]; } calc(k2,p);
for(i=0;i<4;i++){ p[i] = w[i] + 0.5 * dt * k2[i]; } calc(k3,p);
for(i=0;i<4;i++){ p[i] = w[i] + dt * k3[i]; } calc(k4,p);
for(i=0;i<4;i++){ w[i] = w[i] + dt * ( k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i] ) / 6.0; }
}
printf("%f %f %f %f",w[0],w[1],w[2],w[3]);
return 0;
}