/*    ATK IV - Numerical Programming (2006)      */
/*     Excercise 6.4 -  Roots of Bessel J0       */

#include<stdio.h>
#include<math.h>
#define NRANSI
#include<nr.h>
#define nbmax 20

/* You can use the fact that dJ0(x)/dx = -J1(x) */

float f(float x){return bessj0(x);}
void fd(float x, float *y, float *dy){
  *y=f(x);
  *dy=-bessj1(x);
}

int main(void){
  int i,nb=nbmax;
  float xacc = 1.0e-6;
  float x1=0.0;
  float x2=20.0;
  float xb1[nbmax+1],xb2[nbmax+1];
  float x;

  zbrak(f,x1,x2,nbmax,xb1,xb2,&nb);
  for(i=1;i<=nb;i++){
    x=rtsafe(fd,xb1[i],xb2[i],xacc);
    printf("%d %f\n",i,x);
  }

  return 0;

}