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+2*x*t+2*x^2+x^4)
}

J <- 100
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 lax wendroff

# 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] <- .5*nu*(1+nu)*u[n-1,j-1]+(1-nu^2)*u[n-1,j]-.5*nu*(1-nu)*u[n-1,j+1]
  }
  plot(x[c(1:(J+1))],u[n,c(1:(J+1))], ylim=c(-.5,1.5),type = "l", col="red")
  lines(x[c(1:(J+1))],uex(x[c(1:(J+1))],dt*(n-1)*matrix(1,1,J+1)))
  title(paste("t = ",format((n-1)*dt,digits=3,nsmall=3)))
  Sys.sleep(.05)
}