u0 <- function(x){
  (x>=.2)*(x<=.4)
}

uex <- function(x,t){
  u0(x-t/(1+x^2))
}

a <- function(x,t){
  (1+x^2)/(1+x*t+2*x^2+x^4)
}

J <- 50
dx <- 1/J
dt <- dx
Nf <- J

x=dx*c(0,c(1:(J+Nf)))

u <- matrix(ncol = J+1+Nf, nrow = Nf+1)
u[1,] <- u0(x)

# método upwind
# esto es para a>0

# aproxima la solución en [0,1]
for (n in c(2:(Nf+1))) {
  u[n,1] <- 0
  for (j in c(2:(J+1+Nf-n))) {
    nu <- a(x[j],dt*(n-1))*dt/dx
    u[n,j] <- (1-nu)*u[n-1,j]+nu*u[n-1,j-1]
  }
  plot(x[c(1:(J+1))],u[n,c(1:(J+1))], ylim=c(0,1))
  lines(x[c(1:(J+1))],uex(x[c(1:(J+1))],dt*(n-1)*matrix(1,1,J+1)))
  Sys.sleep(.05)
}