Coding

Copy (carefully) the numerical recipes routine ``rk4'', rename it ``rk4ind'', and extend it to allow derivs to access the arrays rho, r, dr:

• Add five new arguments to the prototype for ``derivs'' to produce

  C:       void derivs(float x, float y[], float dydx[],
                       int k, float rho[], float r[], float dr[], int Nu)
  
  Fortran: subroutine derivs(x,y,dydx,k,rho,r,dr,Nu)
           implicit none
           real x,y(*),dydx(*),rho(Nu+1),r(Nu+1),dr(Nu+1)
           integer k,Nu

• Add six new arguments to the ``rk4ind'' prototype and update its prototype for derivs to produce

  C:       void rk4ind(float y[], float dydx[], int n, float x, float h, float yout[],
                 void (*derivs)(float, float [], float [], int, float[], float[], float[], int),
                 int k, int dk, float rho[], float r[], float dr[],int Nu)
  
  Fortran: subroutine rk4ind(y,dydx,n,x,h,yout,derivs,k,dk,rho,r,dr,Nu)
           implicit none
           real y(n),dydx(n),x,h,yout(n),rho(Nu+1),r(Nu+1),dr(Nu+1)
           integer n,k,dk,Nu
           external derivs

• Within ``rk4ind'', treat the indices k and dk as integer versions of x and h: (1) when rk4ind defines xh=x+hh with hh=h/2, define kh=k+dk/2, (2) whenever rk4ind calls derivs add the five new arguments to the call, (3) make sure that the appropriate value for k (either k, kh or k+dk) is passed to each call to derivs, paralleling how rk4ind treats x.