posfun.R

f1 <- expression(eps/(2-x/eps))
f2 <- expression(eps*log(exp(x/eps)+1))
flist <- list(f1,f2)


dfun <- function(x,e,n=0,eps=0.001) {
    alt_less <- switch(n+1,
                  eval(e),
                  eval(D(e,"x")),
                  eval(D(D(e,"x"),"x")))
    alt_gtr <- switch(n+1,
                      x,
                      1,
                      0)
    ifelse(x>eps,alt_gtr,alt_less)
}
    
xvec <- seq(-0.002,0.002,length=51)
res <- list()
for (i in seq_along(flist)) {
    for (n in 0:2) {
        res <- c(res,list(dfun(xvec,flist[[i]],n=n)))
    }
}
res <- unlist(res)    
dd <- expand.grid(i=seq_along(xvec),n=0:2,f=seq_along(flist))
dd$res <- res
dd$x <- xvec[dd$i]

library(ggplot2); theme_set(theme_bw())
ggplot(dd, aes(x,res,colour=factor(f))) +
    facet_wrap(~n,scale="free") +
    geom_line()

ggsave("posfun.png",width=7,height=4)