I am trying to solve coupled heat conduction and mass diffusion for a spherical particle.
I am solving the PDE’s using Finite difference method , with 4th accuracy for first and second order equations.

I am discretizing the spherical particle in 20 points .
num_points = 20
delta = radius/num_points.

I am creating an array with 20 elements as initial state for solving the equations.
X0 = np.zeros((num_points,))

I am trying the to solve the following equations:

1. For diffusion term.

def points_with_conc(T, dt, t)

          if i != num_points-1:           # For nodes 1 to 19 , from centre to second last node. 
           
     
           
            diff1 = Deff*(16*T[i+1]-T[i+2]-30*T[i]+16*T[i-1]-T[i-2])/12*delta**2         
            diff2 = Deff*(b-1)*(8*T[i+1]-T[i+2]-8*T[i-1]+T[i-2])/(12/delta/r_act)
            diff = (diff1+diff2)/eps                         # Effective diffusion term
             
            dXdt = get_dXdt(X, T[i], T[i+num_points])        # Calling the mass source function    
            reac = (dXdt*rho_s/gamma/44e-3/eps)              # Effective Reaction rate term 
           
        else:                                                 # Node-20 , i.e At the Surface 
           
            cs = solve_cs(T[i], cg, eps)                         # Calling the surface 
           
            diff = (Deff*((16*T[i+1]-cs-30*T[i]+16*T[i-1]-T[i-2])/12*delta**2 + (b-1)*(8*T[i+1]-cs-8*T[i-1]+T[i-2])/(12/delta/r_act)/eps))
            dXdt = get_dXdt(X, T[i], T[i+num_points])                            
            reac = (dXdt*rho_s/gamma/44e-3/eps)

2. For Heat conduction term

if i != num_points-1:                # For nodes 1 to 19 , from centre to second last node. 

#Heat Diffusion Term -    (∂T)/(∂t)=(α*(∂^2 T)/(∂r^2 ) + 2/r α*(∂T)/(∂r) + ∂α/∂r*(∂T)/(∂r))+Q     

           
           
            ai = lambda_i/rho_i/csolid                  # Effective Thermal Diffusivity 
           
            diff1 = ai*(16*T[i+num_points+1]-T[i+num_points+2]-30*T[i+num_points]+16*T[i+num_points-1]-T[i+num_points-2])/12*delta**2    #  diffusion terms as shown in the HD equation above 
            diff2 = ai*(b-1)*(8*T[i+num_points+1]-T[i+num_points+2]-8*T[i+num_points-1]+T[i+num_points-2])/12/delta/r_act
            
            diff = diff1+diff2                                     # Average diffusion term
                         
            dXdt = get_dXdt(X, T[i], T[i+num_points])           # Calling the mass source function 
            reac = dXdt*delta_rh/csolid/gamma*rho_s/rho_i        # Effective reaction rate term
           
        else:                                          # Node-20 , i.e At the Surface 
           
            rho_i = X*rho_carb+(1-X)*rho_calc
            lambda_i = (lambda_calc*lambda_carb)/(V_carb*lambda_carb + V_calc*lambda_calc)
           
            ai = lambda_i/rho_i/csolid
           
            Ts = solve_Ts(T[i+num_points], Tg, lambda_i)       # Calling of surface temp function
           
            diff1 = ai*(16*T[i+num_points+1]-Ts-30*T[i+num_points]+16*T[i+num_points-1]-T[i+num_points-2])/12*delta**2    #  diffusion terms as shown in the HD equation above 
  diff2 = ai*(b-1)*(8*T[i+num_points+1]-Ts-8*T[i+num_points-1]+T[i+num_points-2])/12/delta/r_act

            diff = float(diff1+diff2)
            dXdt = get_dXdt(X, T[i], T[i+num_points])
            reac = float(dXdt*delta_rh/csolid/gamma*rho_s/rho_i)

I am getting the following error, when i execute this code.

diff1 = ai*(16*T[i+num_points+1]-T[i+num_points+2]-30*T[i+num_points]+16*T[i+num_points-1]-T[i+num_points-2])/12*delta**2 # diffusion terms as shown in the HD equation above

IndexError: index 40 is out of bounds for axis 0 with size 40

Need some help solving this.
Thanks in advance.