library(Matrix)
library(mgcv)

dt <- 0.001
N <- 20
dx <- 1/20
r <- dt/dx^2
M <- Matrix(0, nrow = N-1, ncol = N-1, sparse = TRUE)
diag(M) <- 1-2*r
sdiag(M,-1) <- r
sdiag(M,1) <- r
x <- matrix(c(1:(N-1))/N,nrow = N-1, ncol = 1)
xp <- rbind(0,x,1)
source("p1.R")
u0 <- ini(x)
par(mfrow=c(1,1))
plot(xp,rbind(0,u0,0),type = "l",ylim = c(0,.6))
title(paste("t = ",format(n*dt,digits=3,nsmall=3)))
for (n in c(1:100)) {
  u <- M%*%u0
  up <- rbind(0,u,0)
  u0 <- u
  Sys.sleep(.05)
  plot(xp,up,type = "l",xlim=c(0,1),ylim = c(0,.6))
  title(paste("t = ",format(n*dt,digits=3,nsmall=3)))
}