有意思的问题。按理来说这个规模的矩阵生成，就算用for循环遍历往里填数也不应该这么慢，也不知道mma是抽了什么风。
其中一个加速的方法是向量化，思路很简单我就直接上代码了
m = DiagonalMatrix[1 - 0.5 h*k[h*#, h*#] &[Range[n + 1]]];
m[[2 ;; -1, 1]] = -0.5*h*k[h*Range[2, n + 1], 0];
Do[m[[j + 1 ;; -1, j]] = -h*k[h*Range[j + 1, n + 1], h*j], {j, 2,
n + 1}];
m = SparseArray@m;
效果如下
然而这种方法其实并不能解决问题，因为当n=20000的时候这样生成系数矩阵会爆内存（我16G内存的电脑直接卡死机了）
事实上作为一个下三角矩阵，想要求解的话根本没有必要吧完整的矩阵生成出来，直接一行一行的迭代就好了，代码如下：
ans = Module[{x, alist},
x = ConstantArray[0., n + 1];
x[[1]] = b[[1]]/(1 - 0.5 h*k[h, h]);
Do[
alist = -h*k[h*i, h*Range[i - 1]];
alist[[1]] = -0.5*h*k[h*i, h*i];
x[[i]] = (b[[i]] - alist.x[[1 ;; i - 1]])/(
1 - 0.5 h*k[i h, i h]), {i, 2, n + 1}];
x]; // AbsoluteTiming
n=4000时效果如下图，可以看到速度快乐很多，而且内存占用也从O(n^2)降到了O(n)。按照这种做法n=20000时用时也只有3秒多
