[R-sig-dyn-mod] bvp chemical reactor?

I found out that if I define the function this way (use absolute value 
of y[1]) then the result will be correct and both solvers will work.
Could someone try to explain this to me please.

Best Regards,
Radovan

THIS WILL WORK!
-----------------------
fun <- function(x, y, p) {
  dy1 <- y[2]
  dy2 <- Pe*Da*abs(y[1])^n + Pe*y[2]
  return(list(c(dy1,dy2)))
}
--------------------------
On 5/5/2016 2:51 PM, Radovan Omorjan wrote:

Hello,

This one is also one of my standard examples (chemical reactor) of 
two point boundary problem.
I tried the three functions (bvptwp(), bvpcol(), bvpshoot()) with the 
default parameters and only bvpshoot() gave the results.
Many other examples gave me no problem in solving them with all the 
three solvers, but this particular one.

QUESTION: I would be very grateful if someone could take a look at 
this example and give me a hint if there is a way to get the result 
with the other two solvers.
Just to mention that I do not have to much experience with these 
solvers and I do not know how to explain my students why these two do 
not work for this particular example.

Best Regards,
Radovan

---------------------------
require(bvpSolve)

fun <- function(x, y, p) {
  dy1 <- y[2]
  dy2 <- Pe*Da*y[1]^n + Pe*y[2]
  return(list(c(dy1,dy2)))
}

Pe <- 1
Da <- 2
n <- 0.5

# initial and ending conditions
init <- c(y <- 1, dy <- NA)
end <- c(y <- NA, dy <- 0)

# Solve bvp
#
# This will not work
#
# sol  <- bvptwp(yini = init, x = seq(0,1,len=50),
#               func = fun, yend = end)
# Error in bvpsolver(1, yini, x, func, yend, parms, order, ynames, 
xguess,  :
#                     The Expected No. Of mesh points Exceeds Storage 
Specifications.
# sol  <- bvpcol(yini = init, x = seq(0,1,len=50),
#               func = fun, yend = end)
# Error in bvpsolver(2, yini, x, func, yend, parms, order, ynames, 
xguess,  :
#                     The collocation matrix is singular for the 
final continuation problem
# This one will work
#
sol  <- bvpshoot(yini = init, x = seq(0,1,len=50),
               func = fun, yend = end, guess=-1)

x <- sol[,1]
y <- sol[,2]
par(mar=c(4,4,1.5,1.5),mex=.8,mgp=c(2,.5,0),tcl=0.3)
plot(x,y)