why this function get wrong result at some larger variable value

I'm using gnu scientific library(gsl) to define parabolic cylinder U function. Here is the definition of this U function. But there is a little difference between my definition and the Mathworld's, cause I want to use mine to use the Gauss-Hermite integration method to calculate the infinity interval integration. So I dismiss the exponential part. Here is my code:

#include <gsl/gsl_errno.h>
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_sf.h>
#include <math.h>
#include <string.h>

double ParabolicCylinderU(double a, double x)
{
   double res;
   double tmp1,tmp2;
   tmp1=sqrt(M_PI)/(gsl_sf_gamma(3./4.+a/2.)*pow(2,a/2.+1./4.))*gsl_sf_hyperg_1F1(a/2.+1./4.,1./2.,pow(x,2)/2.);
   tmp2=sqrt(M_PI)/(gsl_sf_gamma(1./4.+a/2.)*pow(2,a/2.-1./4.))*x*gsl_sf_hyperg_1F1(a/2.+3./4.,3./2.,pow(x,2)/2.);
   res=tmp1-tmp2;
   return res;
}

I use scipy.special.pbdv to examine whether my definition is correct. For the smaller x value, the result is agreed with the pbdv, but when the value gets larger, the result is more weird. How can I solve this problem. There is some part of output:

5.3599999999999994        10.720004342260813
5.4000000000000004        10.801040023779478     
5.4399999999999995        10.879249134730191     
5.4800000000000004        10.961176175487582     
5.5199999999999996        11.036212437780495     
5.5600000000000005        11.115912986845069     
5.5999999999999996        11.189683877125612     
5.6400000000000006        11.289942688341265     
5.6799999999999997        11.385711218186625     
5.7200000000000006        11.379394963160701     
5.7599999999999998        11.532254417763568     
5.8000000000000007        11.575985086141165     
5.8399999999999999        11.881533311150061     
5.8800000000000008        11.896642911301599     
5.9199999999999999        11.639708082791794     
5.9600000000000009        11.550981603033073     
6.0000000000000000        10.455916671990396     
6.0399999999999991        9.3887694230651952 
...
9.6799999999999997        5.1646711481042656E+024
9.7199999999999989       -1.5183962001815768E+025
9.7600000000000016       -3.5383625277571119E+025
9.8000000000000007        4.2306292738245936E+025
9.8399999999999999       -1.3554993237535422E+026
9.8799999999999990       -3.9599339761206319E+026
9.9200000000000017        4.4052629528738374E+026
9.9600000000000009        3.7902904410228328E+027
10.000000000000000       -1.6696086530684606E+021

the first column is 'x' value, second column is parabolic function value, it happens to be wrong at about 5.83999999, because if we input this ParablicCylinderU(-0.999999-1./2.,sqrt(2)*x), this is nearly a straight line. The following is my examination python code.

#!/usr/bin/env python                                                                                                                                                                          

from math import *
import numpy as np
from scipy import special
import matplotlib as mpl
mpl.use('Agg')
import matplotlib.pyplot as plt

nz= 0.999998779431

x=np.linspace(0,10,250,endpoint=True)
d,dv=special.pbdv(nz,sqrt(2.)*(x))

plt.figure()
d=d*np.exp(np.power(x,2)/2)*np.power(2,nz/2.)
plt.plot(x,d)
plt.plot(x,np.zeros(x.shape[0]))
plt.savefig("scipy_parabolic.png")

Solution

You can get pbdv source code from specfun scipy code Look for SUBROUTINE PBDV, if you not understand fortran, you can convert it to C via f2c