Solving a system of partial differential equations (PDE's) in R

Question

I need to solve the following system of PDE's that contains diffusion terms in R:

I use the R package ReacTran to solve the above system. This is my R code:

library(ReacTran)
library(deSolve)
library(rgl)

# Parameters
r <- 1
epsilon <- 1
delta <- 0.1
a_P <- 5
b_P <- 0.1
a_I <- 0.1
b_I <- 0.01
D_C <- 0.00001
D_I <- 0.001

# Gaussian membership function
gaussmf <- function(x, mu, sigma) { 
  return(exp(-((x - mu) ^ 2) / (2 * (sigma ^ 2))))
}

# Discretization of spatial variable
xgrid <- setup.grid.1D(x.up = 0, x.down = 2, N = 200)

# System of PDE's
system<-function(t, state, parameters) {
  with(as.list(c(state, parameters)),{

    dC <- ((r * C) / (epsilon * C + 1)) * (P / (I + 1)) - delta * C + tran.1D(C = C, D = D_C, dx = xgrid)$dC
    dP <- (a_P * C) - (b_P * P)
    dI <- (a_I * C) - (b_I * I) + tran.1D(C = I, D = D_I, dx = xgrid)$dC

    list(c(dC, dP, dI))})
}

# Setting parameters
parameters <- list(r = r, epsilon = epsilon, delte = delta, a_P = a_P, b_P = b_P, a_I = a_I, b_I = b_I)

# Setting initial values

C <- gaussmf(xgrid$x.mid, mu = 1, sigma = 0.02)
P <- rep(0, times =  xgrid$N)
I <- rep(0, times =  xgrid$N)

state <- c(C = C, P = P, I = I)

# Setting time steps
times <- seq(0, 100, by = 1)

# Solving the system
out <- ode.1D(y = state, times = times, func = system, parms = parameters, dimens = xgrid$N)

# Plotting the result
persp3d(x=times, y=xgrid$x.mid, z=out[,2:(xgrid$N+1)],col = "skyblue",xlab = "Time", ylab = "Location", zlab = "C")

I have already solved this system via MATLABs pdepe function and I have get desired result. But by running the above code in R I do not get the right result (like MATLAB). It seems the diffusion term is not working. What is wrong with it?

Edit:

In fact I should get following curves for different values of the parameter a_I:

For a_I = 3: