[R-sig-dyn-mod] bvpshoot - is this a bug?

[Radovan Omorjan]

Hello,

I tried the example from the help of bvpshoot. Just slightly changed the 
code and got two different results.
I think this must be a bug here, or I am doing something wrong.

Please take a look at this code. It will give you two different results 
for the same problem

Regards,
Radovan

Code

## 
=============================================================================
## Example 2 - initial condition is a function of the unknown x
## tubular reactor with axial dispersion
## y’’=Pe(y’+Ry^n) Pe=1,R=2,n=2
## on the interval [0,1] and with initial conditions:
## y’(0)=Pe(y(0)-1),
## y’(1)=0
##
## dy=y2
## dy2=Pe(dy-Ry^n)
## 
=============================================================================
require(bvpSolve)

###
### WRONG ####
###

Reactor<-function(x, y, parms){
   with(as.list(c(parms)),
        {
          list(c(y[2],
                 Pe * (y[2]+R*(y[1]^n))))
        })
}

par <- c(Pe = 1, R = 2, n = 2) ####try to put n=0.5#######

xx <- seq(0, 1, by = 0.01)
# 1. yini is a function here
yini1 <- function (x, parms) {
                     with(as.list(c(parms)),c(x, Pe*(x-1)))
}
system.time(
   sol1 <- bvpshoot(func = Reactor, yend = c(y = NA, dy = 0),
                   yini = yini1, x = xx, extra = 1, parms=par)
)
diagnostics(sol1)
plot(sol1, which = "y", main = "Reactor", type = "l", lwd = 2)
attributes(sol1)$roots
head(sol1)
############################
###
### CORRECT ####
####

Reactor1<-function(x, y, parms)
{
   list(c(y[2],
          Pe * (y[2]+R*(y[1]^n))))
}

Pe <-  1
R <-  2
n <-  2 ##### Try to put n <- 0.5 #########

xx <- seq(0, 1, by = 0.01)
# 1. yini is a function here
yini <- function (x, parms) c(x, Pe*(x-1))

system.time(
   sol <- bvpshoot(func = Reactor1, yend = c(y = NA, dy = 0),
                   yini = yini, x = xx, extra = 1)
)
diagnostics(sol)
plot(sol, which = "y", main = "Reactor", type = "l", lwd = 2)
attributes(sol)$roots
head(sol)
##############################################