Actually it just the parameterization that is causing trouble near k=0
let u = (x-z)/a then the problematic part of your function is (1- k*u)^(1/k) take the log to get log(1-k*u)/k = -(k*u +k^2*u^2/2 + ...)/k = -(u +k u^2/2 + ..) so your function is exp(-u - ku^2/2 - ...) and use that for k<eps for some small eps. You can add the right number of terms and fix the size of eps so that the function behaves well. ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.