Eu gostaria de saber como escrever uma matriz tridiagonal NxN. Eu tenho três vetores que serão atualizados em um loop.

Eu estou trabalhando com algo semelhante a isso:
[![Example][1]][1]

Nota: Eu apenas gostaria de saber sobre zeros e que parâmetros do Python eu poderia usar para ajustar isso. Bem, eu tenho dois códigos aqui, o primeiro eu escrevi em Fortran e está funcionando bem. O segundo é o que eu tentei escrever em Python.

Fortran:

    do i=2,n-1
      do j=2,n-1
        if (i.eq.j) then
    
          D(i,j+1)=-u_med(i+1)/(delta_r(i)*delta_r(i+1))
    
          t1 =u_med(i+1)/(delta_r(i)*delta_r(i))
    
          t2 = u_med(i)/(delta_r(i)*delta_r(i))
    
          D(i,j)= t1 + t2 + V(i)
    
          D(i+1,j)=-u_med(i+1)/(delta_r(i)*delta_r(i+1))
    
        end if
    
      end do
    end do

Python:

    for i in range(2,n):
      for j in range(2,n):
        if i == j:
          D[i][j+1] = - u_med[i+1]/(delta_r[i]*delta_r[i+1])
          t1 = u_med[i+1]/(delta_r[i]*delta_r[i])
          t2 = u_med[i]/(delta_r[i]*delta_r[i])
          D[i][j]= t1 + t2 + V[i]
          D[i+1][j]= - u_med[i+1]/(delta_r[i]*delta_r[i+1])
    
    t1 = u_med[2]/(delta_r[1]*delta_r[1])
    t2 = 0
    D[1][1]= t1 + t2 + V[1]
    D[1][2]= - u_med[2]/(delta_r[1]*delta_r[2])
    D[2,1]= - u_med[2]/(delta_r[2]*delta_r[1])
    t1 = 0
    t2 = u_med[n]/(delta_r[n]*delta_r[n])
    
    D[n][n]= t1 + t2 + V[n]

Que dá o erro:

> D[i][j+1] = - u_med[i+1]/(delta_r[i]*delta_r[i+1]) ValueError: setting
> an array element with a sequence.

Definições de `u_med`, `delta_r`, `r`, `n` e `D` em Python:

    u_med = np.zeros((n,2))
    for i in range(2,n):
      tta1 = r[i]*u[i]
      tta2 = r[i-1]*u[i-1]
      u_med[i] = 0.5*(tta1 + tta2)/(r[i] - r[i-1])
    
    u_med[1] = u_med[2]
    
    delta_r = np.zeros((n-1,2))
    for i in range(2,n-1):
      ft1 = r[i+1]*r[i+1]
      ft2 = r[i-1]*r[i-1]
      ft3 = 2*r[i]*(r[i+1] - r[i-1])
      delta_r[i] = math.sqrt(0.125*abs(ft1 - ft2 + ft3))
    
    ri=0
    n1 = 51
    
    r1 = ri
    r2 = 250
    hr1 = (r2-r1)/(n1-1)
    
    r = np.zeros((n1,1))
    for i in range(n1):
      r[i] = r1 + i*hr1
    
    u = np.zeros((n+1,1))
    for i in range(1,n+1):
      i = 1
    
    D = npm.zeros((n,n))

  [1]: https://i.sstatic.net/t5t7d.png

Eu suspeito que o erro esteja em:

    D = npm.zeros((n,n))

Mas como eu poderia reescrever isso?