Not sure what to make of this

I have the following two functions.  The first is a solution to Shrodinger for the m*k^2*x^4 potential.  The second tries to find an eigenvalue to the first given two initial guess's.  It is not elegant but it seems to work.  I keep getting the same answer for the eigenvalues accurate to 32 number of digits.  The accurate part I don't know for sure.  All I know is I get the same 32 digits all the time.  How accurate is my answer?  The second function prints out the answer two 32 decimal places and it gives that same exact answer for different parameters and number of iterations as long as it is high enough.  Can anyone tell me if this is right, the method seems good but when I change a digit in the answer is has no effect on my 'tests' for correctness.  I suspect the first 20 digits or so may be right, but how to tell is a problem.  Any ideas

ratprint:false;

ratepsilon:10^-100;

terms:1600;

fpp:500;

f(x,y,terms):=
block
([c,p,n,hbar,k,m],
    c[-4]:0,c[-3]:0,c[-2]:0,c[-1]:0,c[0]:1,c[1]:0,hbar:1,m:2,k:2/3,
    for n: 2 thru terms do (c[n]:(k^2 * m^2 * c[n-6] - 2 * m * c[n-2] * y)/(n * (n - 1) * hbar^2)),
    p:c[0],
    for n: 1 thru terms do (p:p + c[n]*x^n),
    p
)$

fsolve(py1, py2, xg, iterations):=
block
(
    fpprec:500,i:0, ii:0, j:0, px:xg, pyb:bfloat(py1), pyt:bfloat(py2),
    while ii < iterations do 
        block
        (
            ft:f(px,pyt,terms),
            st:block(dfx:(f(px+1/100,pyt,terms)-ft)*100, 
                dfy:(f(px,pyt+(pyt-pyb)/1000000,terms)-ft)/((pyt-pyb)/1000000), -dfx/dfy),
            fb:f(px,pyb,terms),
            sb:block(dfx:(f(px+1/100,pyb,terms)-fb)*100, 
                dfy:(f(px,pyb+(pyt-pyb)/1000000,terms)-fb)/((pyt-pyb)/1000000), -dfx/dfy),
            if st=sb then go(endloop),
            p2:(pyb*st-pyt*sb)/(st-sb),
            pyt:(pyt+p2*999)/1000,
            pyb:(pyb+p2*999)/1000,
            ii:ii+1,
            i:ii,
            go (again),
            endloop,
            i:ii,
            ii:10000,
            again,
            9
        ),
        fpprec:32,
        /* return a list of num iterations and high and low value */
        [i,bfloat(pyt),bfloat(pyb)]
)$

With this 

fsolve(8, 9, 6, 15);

->  [15,8.0342385473705851930757373813198b0,8.0342385473705851930757373813198b0]

Thanks,

Rich

PS You can reduce the terms to 1200 with no effect on answer.  You can change 6 to 7 with no effect at 1600 terms.