Numerical solution of the advection equation with Crank–Nicolson finite difference method

I need to implement a numerical scheme for the solution of the one-dimensional advection equation
$$\frac{partial u}{partial t} + C(x, t) frac{partial u}{partial x} = 0 \\$$
$$ \ C(x,t) = frac{picos(2 pi t) + 3.5}{2pi – pisin(pi x)} \\$$
$$\ 0 leq x leq 1, quad 0 leq t leq 1\\$$

with Crank–Nicolson finite difference method
$$\ frac{u^{j + 1}_{k} – u^{j}_{k}}{tau} + frac{C}{4h}(u^{j + 1}_{k + 1} – u^{j + 1}_{k – 1} + u^{j}_{k + 1} – u^{j}_{k – 1}) = 0\\$$
and boundary conditions
$$\ phi(x) = u(x,0), quad psi_{0}(t) = u(0,t), quad psi_{1}(t) = u(1, t) \\$$

it’s known that the exact solution has the following form
$$\ U(x,t) = cos(pi x) – frac{1}{2}sin(2 pi t) + 2pi x – 3.5t\\$$

The numerical solution on each layer in time is found as the solution of a system of linear equations with a tridiagonal matrix
$$\ Au = b \\$$
$$\ A =
begin{pmatrix}
1 & frac{Ctau}{4h} & 0 & …\
-frac{Ctau}{4h} & 1 & frac{Ctau}{4h} & … \
0 & -frac{Ctau}{4h} & 1 & … \
… & … & … & …
end{pmatrix}\\$$

I’ve write a simple python script, but the result doesn’t match the theory:

I’m sure I’m wrong somewhere, but I don’t see where.
Thanks!