Monte carlo method for uncertainty

Question

enter image description here

enter image description here

Hi, currently I'm working on the question background is stated in the figures. and I need to find the expanded uncertainty at a confidence level of 95% around U and N separately for groups A and R using the Monte Carlo Method.

I tried the code below, and the probability density function of variables is given (as coded in the for loop).

Assumption made: -initial U and N for adopters and refusers -The number of iteration -the t = np.linspace(0, 5,5), which relates to rate of change

the results i get were weird , the UA, UR, NA, NR are getting extremely small, which i don't understand why. Did i make any mistake in the code?

import numpy as np
from scipy.integrate import odeint

# Iteration MCM
n = 1000 

# adopter
UA_initial=200 
NA_initial=20000  

# refusers
UR_initial=20   
NR_initial=2000   

UA_result = []
UR_result = []

NA_result = []
NR_result = []



# Define the ODE system
def system(y, t, w, f, g, l, b):
    U,N=y
    dU_dt = w * f * U * N - g * U - l * U 
    dN_dt = -w * f * U * N + b * N
    return dU_dt, dN_dt


for i in range(n):
    w_A_pdf = np.random.uniform(0,1)
    f_A_pdf  = np.random.beta(5,2)
    g_A_pdf  = np.random.poisson(5)
    l_A_pdf  = np.random.poisson(3)
    b_A_pdf  = np.random.normal(1,2)

    w_R_pdf  = np.random.uniform(0,1)
    f_R_pdf  = np.random.beta(2,5)
    g_R_pdf  = np.random.poisson(3)
    l_R_pdf  = np.random.poisson(1)
    b_R_pdf  = np.random.normal(0,4)  

    t = np.linspace(0, 5,5) 
    y0_A = [UA_initial, NA_initial]
    y0_R = [UR_initial, NR_initial]

    # Solve the ODEs
    solution_A = odeint(system, y0_A, t, args=(w_A_pdf , f_A_pdf , g_A_pdf , l_A_pdf , b_A_pdf ))  
    solution_R = odeint(system, y0_R, t, args=(w_R_pdf , f_R_pdf , g_R_pdf , l_R_pdf , b_R_pdf ))

    # Store result
    UA_result.append(solution_A[:,0])  
    NA_result.append(solution_A[:,1])  
    UR_result.append(solution_R[:,0])  
    NR_result.append(solution_R[:,1])