偏微分方程式 (Subscript[u, t]) = (Subscript[u, xx])を初期条件 u(0,x)=100(0<x<π),境界条件 u(t,0)=u(t,π)=0 のもとでの解法

Clear[f]; f[x_] := 100 /; 0 < x < π; 
    f[x_] := 0 /; x = 0 = π; 
    Plot[f[x], {x, 0, π}, PlotStyle -> Hue[0.0], 
      AxesLabel -> None, Ticks -> {{0, π}, {100}}]

Clear[v]; v[n_] := Table[N[f[π (k/n)]], {k, 1, n - 1}]; 
    v[4]

h = (1/m), k = (1/n)とする．このとき，行列 A=A (m, n)を次のように定義する

Clear[o, z]; o[x_] = 0; z[n_] := Array[o, n - 4]; Clear[k]; 
    k[m_, n_] := (n^2)/m; Clear[r]; 
    r[m_, n_, i_] := 
      Insert[Insert[Insert[z[n], k[m, n], i - 1], 1 - 2  k[m, n], i], 
        k[m, n], i + 1]; Clear[rf]; 
    rf[m_, n_] := 
      Insert[Insert[Insert[z[n], 1 - 2  k[m, n], 1], k[m, n], 2], 0, 3]; 
    Clear[rl]; 
    rl[m_, n_] := 
      Insert[Insert[Insert[z[n], 0, n - 3], k[m, n], n - 2], 1 - 2  k[m, n], 
        n - 1]; Clear[M]; 
    M[m_, n_] := {rf[m, n], Table[r[m, n, i], {i, 2, n - 2}], rl[m, n]}; 
    Partition[Flatten[M[10, 6]], 10]; 
    A[m_, n_] := Partition[Flatten[M[m, n]], n - 1] // N

U (m, n, s) の定義

Clear[L]; L[m_, n_, s_] := MatrixPower[A[m, n], s]; Clear[U]; 
    U[m_, n_, s_] := L[m, n, s] . v[n]; U[40, 4, 0]

時間    t = (s/m)における u (t, x)の近似グラフ
　　　　　　　　　u[m, n, s]
             の定義

Clear[u]; 
    u[m_, n_, s_] := 
      Insert[Insert[U[m, n - 1, s], 0, n - 1], 0, 1]

ここで    n,mを ((n^2)/m) < (1/2)を満たすように選んで ,\(時間を\ t\)=(s/m) ( s = 0,1,2, ....)として 
                                                u[m, n, s]
               の挙動を見る．

Do[ListPlot[u[100, 6, s], PlotJoined -> True, PlotStyle -> Hue[0.0], 
        PlotRange -> {0, 105}, AxesLabel -> {x, u}, 
        Ticks -> {{1}, {100}}], {s, 0, 30}]