restart: Digits:=16; N:=8; wt:=exp(-50*(x^4/4-x^2/2)); g[0]:=evalf(int(wt,x=-infinity..infinity),16); beta[0]:=0; Q[0]:=1; mu[0] := g[0]; Q[1]:=x; for i from 2 to 15 do g[i-1]:=evalf(int(wt*Q[i-1]^2,x=-infinity..infinity),16); beta[i-1]:=g[i-1]/g[i-2]; mu[i]:=evalf(int(wt*x^i,x=-infinity..infinity)); Q[i]:=simplify(x*Q[i-1] - beta[i-1]*Q[i-2]): end do; data:=matrix(15,1); for i from 0 to 14 do data[i+1,1]:=beta[i]; end do; writedata("\mydata.dat",data,float);