for ( k in 1:maxk ){
for ( i in 2:( length(U[,1])-1 ) ){
U[i,k+1] <- U[i,k] + 0.5*r*( U[i-1,k] + U[i+1,k] - 2*U[i,k]) + dt*U[i,k]*(1-U[i,k])^2
  for ( i in 1:1){ #ghost sites
    U[i,k+1] <- U[i+2,k+1]
     for ( i in length(U[,1]-1):length(U[,1]-1)){ #ghost sites
       U[i,k+1] <- U[i-2,k+1]
     }
    }
   }
u_new <- U[,k+1]
}

U <- matrix ( , nrow = length ( u ) , ncol = length ( time ))
U[,1] <- u #boundary condition, defined earlier
r<-2*(dt/(dx)^2)