# チューリングパターン　反応拡散系

```a <- 2.8 * 10^(-4)
b <- 5 * 10^(-3)
tau <- 0.1
k <- -0.005

size <- 100
dx <- 2/size
T <- 10.0
dt <- 0.9 * dx^2/2
n <- floor(T/dt)

U <- matrix(runif(size^2),size,size)
V <- matrix(runif(size^2),size,size)

ctr <- 2:(size-1)
plus <- 3:size
minus <- 1:(size-2)

U.hx <- list()
V.hx <- list()
U.hx[[1]] <- U
V.hx[[1]] <- V
for(i in 1:n){
deltaU <- (U[minus,ctr] + U[plus,ctr] + U[ctr,minus] + U[ctr,plus] - 4 * U[ctr,ctr])/(dx^2)
deltaV <- (V[minus,ctr] + V[plus,ctr] + V[ctr,minus] + V[ctr,plus] - 4 * V[ctr,ctr])/(dx^2)

Uctr <- U[ctr,ctr]
Vctr <- V[ctr,ctr]

new.U <- U
new.V <- V
new.U[ctr,ctr] <- Uctr + dt * (a*deltaU + Uctr - Uctr^3 -Vctr + k)
new.V[ctr,ctr] <- Vctr + dt * (b*deltaV + Uctr - Vctr)/tau

new.U[1,] <- new.U[2,]
new.U[size,] <- new.U[size-1,]
new.U[,1] <- new.U[,2]
new.U[,size] <- new.U[,size-1]
new.V[1,] <- new.V[2,]
new.V[size,] <- new.V[size-1,]
new.V[,1] <- new.V[,2]
new.V[,size] <- new.V[,size-1]
U <- new.U
V <- new.V

#U.hx[[i+1]] <- U
#V.hx[[i+1]] <- V

if(i%%5000 == 1){
image(U,col=topo.colors(100))
}
}
```