PDE imaginary

## =============================================================================
## R-Code to solve example 9.3.4 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 partial differential equation model with complex numbers
library(deSolve)

alf <- 0.5
gam <- 1

Schrodinger <- function(t, u, parms) {
   du <- 1i * tran.1D (C = u, D = 1, dx = xgrid)$dC  +
                1i * gam * abs(u)^2 * u
   list(du)                                   
}
N     <- 300    
xgrid <- setup.grid.1D(-20, 80, N = N)    
x     <- xgrid$x.mid                  

c1  <- 1
c2  <- 0.1

sech    <- function(x) 2/(exp(x) + exp(-x))
soliton <- function (x, c1)
  sqrt(2*alf/gam) * exp(0.5*1i*c1*x) * sech(sqrt(alf)*x) 

yini <- soliton(x, c1) + soliton(x-25, c2)

times <- seq(0, 40, by = 0.1)
print(system.time(
 out <- ode.1D(y = yini, parms = NULL, func = Schrodinger, 
              times = times, dimens = 300, method = "adams")
))

image(abs(out), grid = x, ylab = "x", main = "two solitons", 
  legend = TRUE)