#include <stdio.h>
#include <math.h>

#define NumVars	4

double fext(double t)
{
  if (0.00 <= t && t <= 10.00) {
    return 1.00;
  } else {
    return 0.00;
  }
}

void differential(int n, double x[], double t, double value[])
{
  /*	x[0]:x	x[1]:v	x[2]:xspring	x[3]:xvoigt	*/
  
  static double m = 0.1;
  static double k1 = 2.5;
  static double b1 = 2.0;
  static double k2 = 0.2;
  static double b2 = 0.3;
  
  value[0] = x[1];
  value[1] = (-k1*x[2] + fext(t))/m;
  value[2] = x[1] - (k1/b1 + k1/b2)*x[2] + (k2/b2)*x[3];
  value[3] = (k1*x[2] - k2*x[3])/b2;
}

main()
{
  int		kcnt;
  double	t, dt;
  double	x[NumVars];
  char		*filename = "fourelm.dat";
  FILE		*fd;
  
  if ((fd = fopen(filename,"w")) == NULL) {
    fprintf(stderr, "cannot open %s\n", filename);
    exit(1);
  }
  
  kcnt = 20000;
  t = 0.0;
  dt = 0.001;
  
  /* initial values */
  x[0] = 0.0;
  x[1] = 0.0;
  x[2] = 0.0;
  x[3] = 0.0;
  RungeKutta(fd, NumVars, x, differential, dt, kcnt);
  
  fclose(fd);
}