OutOfMemoryError() when approximating a 3D function

[weymouth]

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?