BVP

## =============================================================================
## R-Code to solve example 11.3 from the book
## K. Soetaert, J. Cash and F. Mazzia, 2012.  
## Solving differential equations in R. UseR, Springer, 248 pp. 
## http://www.springer.com/statistics/computational+statistics/book/978-3-642-28069-6.
## implemented by Karline Soetaert
## =============================================================================

## A boundary value problem
library(bvpSolve)
swirl <- function (t, Y, eps)  {
  with(as.list(Y), 
    list(c((g*f1 - f*g1)/eps,
          (-f*f3 - g*g1)/eps))
      )
}

eps  <- 0.001
x    <- seq(from = 0, to = 1, length = 100)
yini <- c(g = -1, g1 = NA, f = 0, f1 = 0, f2 = NA, f3 = NA)
yend <- c(1, NA, 0, 0, NA, NA)
Soltwp <- bvptwp(x = x, func = swirl, order = c(2, 4),
                 par = eps, yini = yini, yend = yend)

pairs(Soltwp, main = "swirling flow III, eps=0.001", 
  pch = 16, cex = 0.5)

diagnostics(Soltwp)