..................の間がプログラム本体です
program 16
..................
/* 
一次元調和振動子の運動方程式をEuler法で計算する
（エネルギーが発散するので、計算法としてはよくない(^^)）
*/
#include<stdio.h>

void euler(double *, double *, double);

main()
{
 double delta_t,q,p;
 int i,step;

/* 時間刻み、ステップ数、初期の座標と運動量を入力 */
 scanf("%lf%d%lf%lf",&delta_t,&step,&q,&p);

/* 1ステップ分を計算する関数をステップ数だけ繰り返して呼ぶ */
 for(i=1; i<step; ++i){
  euler(&q, &p,delta_t);

/* 各ステップごとに、時刻、座標、運動量、エネルギーを出力 */
  printf("%f  %f   %f   %f\n",i*delta_t,q,p,(p*p+q*q)/2);
 }
}

/* Euler法1ステップ分の計算をする関数 */
void euler(double *q, double *p, double delta_t)
{
 double t;            
 t = *q; 
 *q += *p *delta_t;   
 *p += -t *delta_t;
}
..................
program 17
..................
/* 
16と同じ問題を四次のRunge-Kuttaで解くプログラム
*/
#include<stdio.h>

void rungekutta(double *, double *, double);
void onestep(double, double, double *, double *, double);

main()
{
 double delta_t,q,p;
 int i,step;

/* 時間刻み、ステップ数、初期の座標と運動量を入力 */
 scanf("%lf%d%lf%lf",&delta_t,&step,&q,&p);

/* 1ステップ分を計算する関数をステップ数だけ繰り返して呼ぶ */
 for(i=1; i<step; ++i){
  rungekutta(&q,&p,delta_t);

/* 各ステップごとに、時刻、座標、運動量、エネルギーを出力 */
  printf("%f  %f   %f   %f\n",i*delta_t,q,p,(p*p+q*q)/2);
 }
}

/* 1ステップ分の計算をする関数 */
void rungekutta(double *q, double *p, double delta_t)
{ 
  double k1q,k1p,k2q,k2p,k3q,k3p,k4q,k4p;
  onestep(*q, *p, &k1q, &k1p,delta_t);
  onestep(*q+k1q/2, *p+k1p/2, &k2q, &k2p, delta_t);
  onestep(*q+k2q/2, *p+k2p/2, &k3q, &k3p, delta_t);
  onestep(*q+k3q, *p+k3p, &k4q, &k4p, delta_t);
  *q += k1q/6 + k2q/3 + k3q/3 + k4q/6;
  *p += k1p/6 + k2p/3 + k3p/3 + k4p/6;
  }

/* rungekuttaの一段分 */
void onestep(double q, double p, double *qq, double *pp, double delta_t)
{
 *qq = p*delta_t;
 *pp = -q*delta_t;
}