I am trying to obtain a numerical solution to a system of diffusion-advection equations using python. At the moment I am getting a “only size-1 arrays can be converted to python scalars” error on the line just past the double for loop beginning with g[i,j+1]. I am not sure why because all of my arrays are 2D. I am also not sure if my system requires a more advanced numerical method to obtain a good numerical solution due to the source terms and large numbers. Also s depends on g. The equations along with the initial and boundary conditions are in the link below. All parameters are in units of cm, g, and s. The variables g and s stand for gas and stars respectively. I am trying to see how an initial mass distribution for the quantities g and s evolve over some time t which does not have to be a specific time. I greatly appreciate any help on how to proceed.

System of equations with initial and boundary conditions and values of the parameters

import numpy as np
import matplotlib.pyplot as plt

# Time dependence
Nt = 1000          # number of timesteps
dt = 1e7

# Space dependence
Nr = 901           # number of nodes
rmin, rmax = 0.1, 1.0
dr = ( rmax - rmin ) / (Nr - 1)
r= np.linspace(rmin, rmax, Nr)

#Values of parameters
r0 = 1.56e18 #Initial radius
v = 2.08e8  #Velocity
nu=1e19 #Viscosity
m=1e34 #Mass distribution of gas which is also set to equal the mass distribution of stars

#Initializing array to zero
g = np.zeros((Nr,Nt))   
s = np.zeros((Nr,Nt))

# Dependent variables with Gaussian initial conditions
mu = .5*(rmax - rmin)
sigma = 0.1*(rmax-rmin)

g[i,0] = np.exp(-(r[i]-mu)**2/(2*sigma**2)) / np.sqrt(2*np.pi*sigma**2)
s[i,0] = np.exp(-(r[i]-mu)**2/(2*sigma**2)) / np.sqrt(2*np.pi*sigma**2)

for j in range(0,Nt-1):   
    for i in range(1,Nx-1):  
        
        g[i,j+1] = g[i,j] + 3*nu*(dt/dr**2)*(g[i+1,j] - 2*g[i,j] + g[i-1,j]) + 
        (9/2)*nu*(1/r[i])*( (g[i+1,j] - g[i-1,j])/(2*dr) )*dt   -  
        ((v/r[i])*( g + r[i]*( (g[i+1,j] - g[i,j])/(dr) )))*dt
        
        s[i,j+1] = s[i,j] + (v/r[i])*( g + (r[i]* ((g[i,j] - g[i+1,j])/(dr))))*dt 

I first started with just the advection equation and now included diffusion. I am expecting to see diffusion simultaneously with advection for the gas and stars moving towards the left.