This worksheet computes a finite difference approximation to the solution of the one-dimensional wave equation.

Clear[u]; Clear[A]; Clear[B]; Clear[a]; Clear[b]

Start by clearing variables and setting initial conditions.

f[x_] := 
      10^7*Which[x < (3/8), 0, x > (5/8), 0, True, 
          Exp[(1/((x - 0.5)^2 - (1/16)))] ]

g[x_] := 0

Plot[f[x], {x, 0, 1}]

Now fix the boundary conditions.

Do[u[m, 0] := 0, {m, -1, 100}]

Do[u[m, 40] := 0, {m, -1, 100}]

Set the size of the parameter s.

Choose the size of the time step.

tau = (1.0/50)

xi = (1.0/40)

s = tau/xi

For[j = 1, j < 40, j++, 
      u[-1, 
          j] = -tau*g[j*xi] + (1 - s)*
            f[j*xi] + (s/2)*(f[(j + 1)*xi] + 
                f[(j - 1)*xi])]

For[j = 1, j < 40, j++, u[0, j] = f[j*xi]]

Do[Do[
        u[n + 1, j] = 
          2*(1 - s)*u[n, j] + s*(u[n, j + 1] + u[n, j - 1]) - 
            u[n - 1, j], {j, 1, 39}], {n, 0, 100}]

Do[Do[U[n] = Table[u[n, j], {j, 0, 40}]], {n, 0, 100}]

Do[g[j] = 
            ListPlot[U[j], PlotJoined -> True, 
              PlotRange -> {-1.5, 1.5}];, {j, 0, 100}];