library(Matrix)
library(mgcv)
N <- 100
h <- 1/N
epsi <- .001
M <- Matrix(0, nrow = N-1, ncol = N-1, sparse = TRUE) 
diag(M) <- 2*epsi/h^2
sdiag(M,-1) <- -epsi/h^2-1/(2*h)
sdiag(M,1) <- -epsi/h^2+1/(2*h)
b <- c(rep(1,N-1))
u <- solve(M,b)
u0<-rbind(0,u,0)
x<-seq(0,1,1/N)
plot(x,u0,type = "l")