f(x,En,k,m,hbar,k1,k2,max):=do ([p],fpprec:100,
for n:-6 thru max do (c[n]:0),c[0]:k1,c[1]:k2,
for n:2 thru max do
(c[n]:bfloat((c[n-6]*k^2*m^2-2*c[n -2]*En*m)/(n*(n-1)*hbar^2))),
p:c[0] + sum(c[n]*x^n,n,1,max),
if p > 2.5 then return(2.5) else return(float(p)));
showtime:true;
plot3d('f(x,e,2/3,2,1,1,0,120), [x,-3.99,3.99], [e,0,49.99], [plot_format,gnuplot],
[grid,300,300],
[gnuplot_preamble, "set title 'Anharmonic Potential m k^2 x^4 (300 by 300) First 120 Terms';
set xlabel 'Psi^2'; set ylabel 'Energy';set pm3d map; unset surf;"])$
If you change the grid size to 150,150, this works, but I wanted to go higher.