Explicit finite difference method for transport equation

Question

So, here is an equation

$\begin{array}{c} \frac{\partial u}{\partial t}+a \frac{\partial u}{\partial x}=f(x, t), \quad a>0 \\ 0 \leqslant x \leqslant 1, \quad 0 \leqslant t \leqslant 1 \\ u(x, 0)=\varphi(x), \quad x \in[0,1] \\ u(0, t)=g_{1}(t), \quad t \in[0,1] \end{array}$

A scheme for it:

$\begin{array}{c} \frac{U_{j}^{n+1}-U_{j}^{n}}{\tau}+a \frac{U_{j}^{n}-U_{j-1}^{n}}{h}=f_{j}^{n} \\ j=\overline{1, M} ; \quad U_{0}^{n+1}=g_{1}^{n+1} \end{array}$

and there is an analytic solution:

$\begin{aligned} a &=0.23 \\ u(x, t) &=x^{3}-\sin (2 \pi t x) / 2+2.5 t x \end{aligned}$

It is easy to express $U_{j}^{n+1}$:

$U_{j}^{n+1} = U_{j}^{n} - \frac{a\tau}{h}(U_{j}^{n} - U_{j-1}^{n}) + f\tau $

It seems to be very easy to program this stuff, but I am struggling with it.. I just can't understand what is wrong! Solution just weird.. not even close to truth.

Please, check my code and tell me what is wrong

import numpy as np
import matplotlib.pyplot as plt

# wave speed
a = 0.23

# spatial domain
xmin = 0
xmax = 1

n = 10 # num of grid points

# x grid of n points
X, dx = np.linspace(xmin,xmax,n,retstep=True)
T, dt = np.linspace(xmin,xmax,n,retstep=True)

# initial conditions
def initial_u(x):
    return x**3

def analytic_u(x,t):
    return x**3 - np.sin(2*np.pi*t*x)/2 + 2.5*t*x

def f(x,t):
  return 2.5*x - np.pi*x*np.cos(2*np.pi*t*x) - a*(np.pi*t*np.cos(2*np.pi*t*x) + 2.5*t + 3*x**2)

# each value of the U array contains the solution for all x values at each timestep
U = []

# explicit solution
def u(X,t):
    if t == 0: # initial condition
        return initial_u(X)
    uvals = [] # u values for this time step
    for j in range(len(X)):
        if j == 0: # left boundary
            uvals.append(analytic_u(0,T[t-1]))
        elif j == n-1: # right boundary
            uvals.append(analytic_u(1,T[t-1]))
        else:
            uvals.append( U[t-1][j] - (a*dt/dx)*(U[t-1][j]-U[t-1][j-1]) + f(X[j],T[t-1])*dt )
    return uvals

for t in range(10):
    U.append(u(X,t))
    plt.plot(X,U[-1])