/*           ATK IV - Numerical Programming (2006)             */
/*    Excercise 10.3 -  DE's Continued, Step-size control      */

#include<stdio.h>
#include<math.h>
#define NRANSI
#include<nr.h>
#include<nrutil.h>

int kmax,kount;
float *xp,**yp,dxsav;

void f(float x,float y[],float dy[]);

int main(void){
  float a,b,x,h;
  float epsilon = 1.0e-6;

  float y[] = {0.0,0.0,1.0};

  int neq,ok,bad;
  int i,n;

  neq=2;
  a=1;
  b=3;
  n=1000;
  h=(b-a)/(n-1.0);

  dxsav=h;
  kmax=n;
  xp=vector(1,n);
  yp=matrix(1,neq,1,n);

  odeint(y,neq,a,b,epsilon,h,0.,&ok,&bad,f,rkqs);

  for(i=1;i<=kount;i++){
    printf("%f %f %f\n",xp[i],yp[1][i],yp[2][i]);
  }
  free_matrix(yp,1,neq,1,kount);
  free_vector(xp,1,kount);

  return 0;
  
}

void f(float x,float y[],float dy[]){
  dy[1] = y[2];
  dy[2] = (2./x)*y[2]- (1.0+2.0/(x*x))*y[1];
}