OutOfMemoryError() when approximating a 3D function

I'm trying to recreate a 3D Chebyshev surrogate model from an old paper Newman 1987.

using SpecialFunctions,QuadGK
using Base.MathConstants: γ
function bruteN(x,y,z;withS=true)
    ζ(t) = (z*sqrt(1-t^2)+y*t+im*abs(x))*sqrt(1-t^2)
    Ni(t) = imag(expintx(ζ(t))+log(ζ(t))+γ)
    2/π*quadgk(Ni,-1,0,1,atol=eps())[1]+ifelse(withS,S(x,y,z),0)
end
S(x,y,z) = -2*(1-z/(hypot(x,y,z)+abs(x)))
function bruteNsphere(R,θ,α)
	x,ρ = R*sin(θ),R*cos(θ)
	y,z = ρ*sin(α),-ρ*cos(α)
	bruteN(x,y,z,withS=false)
end

using ApproxFun
R,θ,α = Fun(0..2),Fun(CosSpace()),Fun(0..π/2)
f = bruteNsphere(R,θ,α)

As you can see, its an integral over the special function expintx which makes it very expensive to evaluate, but the integral itself is a smooth function of x,y,z or R,θ,α. (Note that I've subtracted off the log singularity at R=0, which is also done in the paper.)

Newman uses a triple Chebyshev product to estimate this function in 1987 (though he said it took him 150 hours to compute the coefficients). However, when I try to run the code above I get an OutOfMemoryError(). Certainly, I've got more memory in 2024 than Newman had in 1987 😀.

I see there is limited support for 3D functions. Is there any work around I can use here?

Start by reproducing the OutOfMemoryError with the provided bruteNsphere and ApproxFun example, focusing on the 3D function construction and the expensive expintx integral. Review the existing limited 3D-function support; done would be a documented, workable approach for approximating this smooth function without exhausting memory.