A = zeros(M+1);
U2 = zeros(M+1,N+1);
A(1,1) = 1;
for i=2:M
A(i,i-1) = -sigma;
A(i,i) = 1+2*sigma;
A(i,i+1) = -sigma;
end
A(M+1,M+1) = 1;
U2(:,1) = linsolve(A,U(:,1));
for i=2,N+1
U2(:,i) = linsolve(A,U2(:,i-1));
end
Note that the code could be made faster by factorizing A once and solving the linear systems for only changing right-hand sides. This is possible since A is constant for all times.
