$ c++ -o run run.cc -lm -lgsl -lgslcblas
$ run
1.1     0.349982        0.183584
1.2     0.243029        0.155227
1.3     0.136783        0.1292
1.4     0.032991        0.105342
1.5     -0.0669617      0.0835034
1.6     -0.161994       0.0635447
1.7     -0.251279       0.0453342
1.8     -0.334207       0.0287488
1.9     -0.41035        0.0136734
2       -0.479431       -1.15521e-16

// calculate Rydberg wavefunction with given quantum defect.
#include <iostream>
#include <cmath>
#include <fstream>
#include <gsl/gsl_sf_coulomb.h>

using namespace std;

int main(void)
{
 double eta,x,L;
 int k=0;
 gsl_sf_result F,Fp,G,Gp,C_L;
 double exp_F,exp_G;
 double F_L,G_L;

 double kk,KE;
 double z1=1.0;
 double z2=-1.0;
 double wavefunc;

 int n=2;
 double r = 1.0;

 int imax = 10;
 double dr = 0.1;

 for( int i=0; i < 10 ; i++){
   r += dr;

   KE = -1.0/(2.0*n*n)+1.0/r;    //kinetic enegry
   kk = sqrt( 2.0*KE );             //impulse
   eta = z1*z2/kk;                   //conversion of variables
   x = kk*r;                            //conversion of variables

   int status1 = gsl_sf_coulomb_wave_FG_e (eta, x, L, k, &F, &Fp, &G, &Gp, &exp_F, &exp_G);

   int status2 = gsl_sf_coulomb_CL_e (L, eta, &C_L);

   F_L = F.val * exp(exp_F);
   G_L = G.val * exp(exp_G);

   wavefunc = C_L.val*F_L;

   double R_nl = gsl_sf_hydrogenicR ( n, (int)trunc(L), 1.0, r);

   cout << r << '\t' << wavefunc << '\t' << R_nl << endl;

 }//i-loop

 return 0;
}