好吧，我仔细查看了代码。我可以给你的第一个建议是尝试理解代码是如何工作的。很容易理解for循环是一个时间循环，即系统在循环的每次迭代中都向前推进一个时间步。时间行进方程使用 theta 方案进行时间离散化（theta=0.5 --> Crank-Nicolson 方案），并使用 Scipy 的稀疏线性求解器求解得到的线性系统。这个线性系统的分辨率给出了向量$u$每个网格点上的离散温度（或您正在计算其演变的实际物理值）。您编写代码的方式，它存储 u 的过去值。所以要么保存（例如，在列表中），要么只使用$u$在您想要的时候随时运行。

顺便说一下更新方法$Nt$错了，以及计算方法$dt$，所以我已经相应地修改了你的代码。

另外，尝试发布直接可用的代码；）我必须导入一些模块，并设置对solver函数的调用以使其工作。

这是调用的最终结果solver(L=1, Nx=30, Nt=400, T=200., theta=1)（所以 30 个元素）：

顺便说一句，模型似乎有问题，因为它没有在空间中收敛。例如，用 3000 分而不是 30 分，我得到下一个数字。而且情况只会随着点数的增加而恶化$N_x$增加了...我不知道您正在建模的是什么，但是应该检查一下。我还看到dx您的代码中没有使用网格间距，这似乎相当令人惊讶。

代码是：

import numpy as np
import matplotlib.pyplot as plt
import scipy.sparse.linalg

def solver(L, Nx, Nt, T, theta=0.5):
    x = np.linspace(-L,L, Nx+1)   # mesh points in space
    dx = (L-(-L))/(Nx-1)
       
    t = np.linspace(0, T, Nt+1)   # mesh points in time
    dt = t[1]-t[0]

    h = (L-(-L))/(Nx+1)

    u   = np.zeros(Nx+1)   # solution array at t[n+1]
    u_n = np.zeros(Nx+1)   # solution at t[n]
    a   = np.zeros(Nx+1) + (1+x**2)
    c   = np.ones(Nx+1) * -1

    Dl = 0.5*theta
    Dr = 0.5*(1-theta)
   
    diagonal = np.zeros(Nx+1)
    lower    = np.zeros(Nx)
    upper    = np.zeros(Nx)
    b        = np.zeros(Nx+1)
 
    diagonal[1:-1] = 1 + Dl*(a[2:] + 2*a[1:-1] + a[:-2])
    lower[:-1] = -Dl*(a[1:-1] + a[:-2])
    upper[1:]  = -Dl*(a[2:] + a[1:-1])
    
    diagonal[0] = 1
    upper[0] = 0
    diagonal[Nx] = 1
    lower[-1] = 0

    A = scipy.sparse.diags(
        diagonals=[diagonal, lower, upper],
        offsets=[0, -1, 1],
        shape=(Nx+1, Nx+1),
        format='csr')
    
    def I(x):
        return(np.sin(np.pi*x))

    u_n[:Nx+1] = I(x[:Nx+1])
      
    # Time loop
    for n in range(0, Nt):
        print('iter {}/{}'.format(n,Nt))
        b[1:-1] = u_n[1:-1] + Dr*(
            (a[2:] + a[1:-1])*(u_n[2:] - u_n[1:-1]) -
            (a[1:-1] + a[0:-2])*(u_n[1:-1] - u_n[:-2])) + dt*theta*c[1:-1] + dt*(1-theta)*c[1:-1]
        
        # Boundary conditions
        b[0]  = 0
        b[-1] = 0
        # Solve
        u[:] = scipy.sparse.linalg.spsolve(A, b)

        u_n,u = u, u_n 
        
        if 0 :# One plot for every selected iteration
          if np.mod(n,100)==0: # every 100 iterations
            plt.figure()
            plt.plot(x,u)
            plt.grid()
            plt.xlabel('x')
            plt.ylabel('u')
            plt.title('Solution at iteration {}'.format(n))
            plt.show()
        else: # one common figure for all the plotted iterations
          if np.mod(n,75)==0: # every 100 iterations
            if n==0: # instantiate the figure
              plt.figure()
              plt.grid()
              plt.xlabel('x')
              plt.ylabel('u')
            # plot the current iteration
            plt.plot(x,u, label='t={:.2f}'.format(t[n]))
            
            
if __name__=='__main__':
  solver(L=1, Nx=3000, Nt=400, T=200., theta=1)
  plt.legend(loc='center left', ncol=2)

编辑： 通过查看您说您改编的原始文件（https://github.com/hplgit/fdm-book/blob/master/src/diffu/diffu1D_vc.py），t 似乎是定义Dl并且Dr没有意义关于扩散算子的离散化。它应该涉及扩散系数和网格间距，而不是theta时间方案的系数。