我的代码有点长，但是它输出了一个像下面这样的图。我想改变x轴的时间值为一些自定义的东西，比如年（2020，2021，2022），并希望它说的每一个地方40000或42000被替换为年。什么是最简单的方法来做到这一点？（注意，在这种情况下，我不能简单地绘制不同的时间值）

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from scipy.optimize import minimize
import math

# Total population, N.
N = 1
N_m = 1
# Initial number of infected and recovered individuals, I0 and R0.
I0, R0 = 0.001, 0
I0_m, R0_m = 0.001, 0
# Everyone else, S0, is susceptible to infection initially.
U0 = N - I0 - R0
U0_m = N_m - I0_m - R0_m
J0 = I0
J0_m = I0_m
Lf0, Ls0 = 0, 0
Lf0_m, Ls0_m = 0, 0
# Contact rate, beta, and mean recovery rate, gamma, (in 1/days).
beta, gamma = 13.21245908, 365/75
beta_m, gamma_m = 13.21245908, 365/75
mu, muTB, sigma, rho = 1/80, 1/6, 1/6, 0.03
mu_m, muTB_m, sigma_m, rho_m = 1/80, 1/6, 1/6, 0.03
u, v, w = 0.88, 0.083, 0.0006
u_m, v_m, w_m = 0.88, 0.083, 0.0006
t = np.linspace(0, 50000, 50000+1)
int_gamma = 365/56

# The SIR model differential equations.
def genpop(y, t, N, beta, gamma, mu, muTB, sigma, rho, u, v, w):
    U, Lf, Ls, I, R, cInc = y
    b = (mu * (U + Lf + Ls + R)) + (muTB * I)
    lamda = beta * I
    clamda = 0.2 * lamda
    dU = b - ((lamda + mu) * U)
    dLf = (lamda*U) + ((clamda)*(Ls + R)) - ((u + v + mu) * Lf)
    dLs = (u * Lf) - ((w + clamda + mu) * Ls)
    dI = w*Ls + v*Lf - ((gamma + muTB + sigma) * I) + (rho * R)
    dR = ((gamma + sigma) * I) - ((rho + clamda + mu) * R)
    cI = w*Ls + v*Lf + (rho * R)
    return dU, dLf, dLs, dI, dR, cI

# Integrate the SIR equations over the time grid, t.
solve = odeint(genpop, (U0, Lf0, Ls0, I0, R0, J0), t, args=(N, beta, gamma, mu, muTB, sigma, rho, u, v, w))
U, Lf, Ls, I, R, cInc = solve.T

# The SIR model differential equations.
def derivint(y, t, N, beta, int_gamma, mu, muTB, sigma, rho, u, v, w):
    U, Lf, Ls, I, R, cInc = y
    int_gamma = (t + 147368.42105277)/30281.182408084 
    b = (mu * (U + Lf + Ls + R)) + (muTB * I)
    lamda = beta * I
    clamda = 0.2 * lamda
    dU = b - ((lamda + mu) * U)
    dLf = (lamda*U) + ((clamda)*(Ls + R)) - ((u + v + mu) * Lf)
    dLs = (u * Lf) - ((w + clamda + mu) * Ls)
    dI = w*Ls + v*Lf - ((int_gamma + muTB + sigma) * I) + (rho * R)
    dR = ((int_gamma + sigma) * I) - ((rho + clamda + mu) * R)
    cI = w*Ls + v*Lf + (rho * R)
    return dU, dLf, dLs, dI, dR, cI

# Integrate the SIR equations over the time grid, t.
solveint = odeint(derivint, (U[-1], Lf[-1], Ls[-1], I[-1], R[-1], J0), t, args=(N, beta, int_gamma, mu, muTB, sigma, rho, u, v, w))
Uint, Lfint, Lsint, Iint, Rint, cIncint = solveint.T

J_diff = cInc[1:] - cInc[:-1]
J_diffint = cIncint[1:] - cIncint[:-1]
#J_diff = np.diff(cInc)
fig = plt.figure(facecolor='w')
ax = fig.add_subplot(111, facecolor='#dddddd', axisbelow=True)
#ax.plot(t, U*100000, 'black', alpha=1, lw=2, label='uninfected')
#ax.plot(t, Lf/100000, 'r', alpha=1, lw=2, label='latent fast')
#ax.plot(t, (Ls+Lf)*100000, 'black', alpha=1, lw=2, label='latent slow')
#ax.plot(t, I*100000, 'green', alpha=1, lw=2, label='infected')
#ax.plot(t, R*100000, 'red', alpha=1, lw=2, label='recovered')
ax.plot(t[1:], J_diff*100000, 'blue', alpha=1, lw=2, label='Baseline')
ax.plot(t[1:]+(40000-1), J_diffint*100000, 'red', alpha=1, lw=2, label='Reduced delay in diagnosis')
#ax.plot(t, cInc, 'red', alpha=1, lw=2, label='Prevalence')
ax.set_xlabel('Time')
ax.set_ylabel('Number')
ax.set_xlim(39990, 50000)
ax.grid(b=True, which='major', c='w', lw=2, ls='-')
legend = ax.legend()
legend.get_frame().set_alpha(0.5)
plt.title("Incidence")
#plt.savefig('filename.png')
plt.show()