wave <- function ( u, alpha, xdelta, tdelta, n )

#*****************************************************************************80
#
## wave solves the 1D wave equation using the forward time centered space method.
#
#  Licensing:
#
#    Copyright 2016 James P. Howard, II
#
#    The computer code and data files on this web page are distributed under
#    https://opensource.org/licenses/BSD-2-Clause, the BSD-2-Clause license.
#
#  Modified:
#
#    30 March 2020
#
#  Author:
#
#    Original R code by James Howard;
#    Modifications by John Burkardt.
#
#  Reference:
#
#    James Howard,
#    Computational Methods for Numerical Analysis with R,
#    CRC Press, 2017
#    ISBN13: 978-1-4987-2363-3.
#
{
  m <- length ( u )
  uarray <- matrix ( u, nrow = 1 )
  newu <- u

  h <- ( alpha * tdelta / xdelta )^2

  oldu <- rep ( 0, m )
  oldu[2:(m-1)] <- u[2:(m-1)] + h * ( u[1:(m-2)] - 2.0 * u[2:(m-1)] + u[3:m] ) / 2.0

  for ( i in 1 : n )
  {
    ustep1 <- 2.0 * u - oldu
    ustep2 <- u[1:(m-2)] - 2.0 * u[2:(m-1)] + u[3:m]
    newu <- ustep1 + h * c ( 0, ustep2, 0 )
    oldu <- u
    u <- newu
    uarray <- rbind ( uarray, u )
  }

  return ( uarray )
}