trd[a_, b_, c_, 
          n_]  :=  (r  =  Join[{b, c}, Table[0, {k, n - 3}], {a}];  
          m = Table[RotateRight[r, k - 1], {k, n}];  
          m[[1, n]] =  m[[n, 1]]  =  0; m);
trd[b_, c_, n_]  :=  trd[c, b, c, n]

IteratePlot[
        k_]  :=  (Do[vv  =  {xs, uu} // Transpose; 
          vv  =  Join[{{0, 0}}, vv, {{1, 0}}];  
          ListPlot[vv, PlotJoined  ->  True,  
            PlotRange  -> {{0, 1}, yh {-1, 1}}, 
            PlotStyle  -> {Thickness[ .005], Hue[hu]}, 
            PlotLabel  ->  lab]; t  =  t + dt;  
          uu  =  M .  uu,  {i, k}])

Heat Equation — Explicit Scheme

Heatex[k_, n_, kt_]  :=  ( 
          hx  = 1/n // N;   μ  =  kt/hx^2; 
          Print["μ =  ", μ]; 
          xs  =  (Range[n - 1]/n  // N); 
          uu  =  f  /@  (Range[n - 1]/n  // N); 
          M = A  =  trd[1 - 2  μ, μ, n - 1]; t = 0; 
          dt  =  kt; hu = Which[μ  <  1,  .6, True,  .1]; 
          lab  =  
            StyleForm[
              StringJoin["Δx = ", 
                hx // ToString, "     Δt = ", 
                kt // ToString], FontSize  ->  20]; 
          IteratePlot[k]);

f[x_]  :=  
      3  Which[x  <   .2, -x,  x  <   .7,  x  -   .4,  True, 
          1 - x]

yh  =  1.8; Heatex[10, 10,  .01]

yh  =  1.8; Heatex[30, 10,  .005]

yh  =  1.8; Heatex[10, 100,  .005]

yh  =  1.8; Heatex[10, 100,  .0005]

yh  =  1.8; Heatex[30, 100,  .0001]

yh  =  1.8; Heatex[100, 100,  .00005]

Heat Equation — Implicit Scheme

Heatim[k_, n_, kt_]  :=  ( 
          hx  = 1/n // N;   μ  =  kt/hx^2; 
          Print["μ =  ", μ]; 
          xs  =  (Range[n - 1]/n  // N); 
          uu  =  f  /@  (Range[n - 1]/n  // N); 
          A  =  trd[1 + 2  μ, -μ, n - 1]; M =  Inverse[A]; 
          t = 0; dt  =  kt; hu =  .85; 
          lab  =  
            StyleForm[
              StringJoin["Δx = ", 
                hx // ToString, "     Δt = ", 
                kt // ToString], FontSize  ->  20]; 
          IteratePlot[k]);

yh  =  1.8; Heatim[30, 10,  .01]

yh  =  1.8; Heatim[10, 10,  .1]

yh  =  1.8; Heatim[30, 100,  .01]

yh  =  1.8; Heatim[10, 100,  .1]

Heat Equation — Crank–Nicholson Scheme

HeatCN[k_, n_, kt_]  :=  ( 
          hx  = 1/n // N;   μ  =  kt/hx^2; 
          Print["μ =  ", μ]; 
          xs  =  (Range[n - 1]/n  // N); 
          uu  =  f  /@  (Range[n - 1]/n  // N); 
          B =  trd[1 + μ, -μ/2, n - 1]; 
A =  
            trd[1 - μ, μ/2, n - 1]; M =  Inverse[B] . A; t = 0; 
          dt  =  kt; hu =  .1; 
          lab  =  
            StyleForm[
              StringJoin["Δx = ", 
                hx // ToString, "     Δt = ", 
                kt // ToString], FontSize  ->  20]; 
          IteratePlot[k]);

yh  =  1.8; HeatCN[30, 10,  .01]

yh  =  1.8; HeatCN[10, 10,  .1]

yh  =  1.8; HeatCN[30, 100,  .01]