I'm trying to solve the following first order hyperbolic PDE problem using method of lines:Hyperbolic Equation: $u_t = -u_x$ with initial condition: $u(0,x) = 0, 0 < x < 1$boundary condition: $u(t,0) = 1, t \ge 1$The solution should be a step function to the right with velocity $1$. I'm using centered finite difference to get an approximation of $u_x$. Following the code in this tutorial, my code becomesimport numpy as npfrom scipy.integrate import odeintimport matplotlib.pyplot as pltfrom scipy.fftpack import diff as psdiffN = 100 #no. of mesh pointsL = 1.0x = np.linspace(0, L, N) #mesh points xi, 0 < xi < 1h = x[1] - x[0]k = -1.0def odefunc(u, t):    ux = np.zeros(x.shape)    u[0] = 1 # boundary condition    for i in range(1,N-1):        ux[i] = float(u[i+1] - u[i-1])/(2*h)     #   ux[i] = float(u[i] - u[i-1])/h     dudt = -ux    return dudtinit = np.zeros(x.shape, np.float) #initial conditiontspan = np.linspace(0.0, 2.0, N)sol = odeint(odefunc, init, tspan, mxstep=5000)for i in range(0, len(tspan), 2):    plt.plot(x, sol[N-1], label='t={0:1.2f}'.format(tspan[i]))plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))plt.xlabel('t')plt.ylabel('u(x,t)')plt.subplots_adjust(top=0.89, right=0.77)plt.savefig('pde.png')from mpl_toolkits.mplot3d import Axes3Dfig = plt.figure()ax = fig.add_subplot(111, projection='3d')SX, ST = np.meshgrid(x, tspan)ax.plot_surface(SX, ST, sol, cmap='jet')ax.set_xlabel('x')ax.set_ylabel('t')ax.set_zlabel('u(x,t)')ax.view_init(elev=15, azim=-100) # adjust view so it is easy to seeplt.savefig('pde-3d.png')However the resulting graph is not a step function as it should be. What might be the problem here? 
Hyperbolic Equation PDE (Python)
pde;finite difference;python;hyperbolic pde
First-order hyperbolic equations model conservation laws; as the alternative name transport equations suggests, they transport information along so-called characteristic curves with a finite speed of propagation.For the simple model equation $u_t + u_x = 0$, you can show (e.g., by separation of variables) that the solution should be of the form$$u(t,x) = u(0,x-t),$$where the missing information (where $x-t<0$) is taken from the boundary conditions. As you correctly write, with the data in your question you'd expect a step profile moving to the right with constant speed $1$:$$u(t,x) = \begin{cases} 1 & x\leq t,\\0& x>t. \end{cases}$$Here you can already see that something should be going wrong if you naively use a finite difference approximation: This solution is not differentiable (it's not even continuous), and yet you try to approximate its derivative using a finite difference quotient. So you need to be a bit more clever.One way is to slightly modify the equation so that it does have a differentiable solution, for example by considering $u_t +u_x = \varepsilon u_{xx}$ for $\varepsilon>0$ very small, and apply a finite difference method to that. It turns out that taking $\varepsilon = \frac{\Delta t}{2}$ (half the discrete time step size) is a good idea; this is called the Lax-Wendroff method:$$U^{n+1}_{j} = U^n_j - \frac{\Delta t}{2\Delta x}(U^n_{j+1}-U^n_{j-1}) + \frac{\Delta t^2}{2\Delta x^2}(U^n_{j+1}-2U^n_j + U^n_{j-1}).$$You can show that this is second-order accurate and stable as long as $\Delta t\leq \Delta x$. (This is discussed in Chapter 10 of Randy LeVeque, Finite Difference Methods for Ordinary and Partial Differential Equations, SIAM 2007, together with other schemes).So if you replace your odefunc routine byt = np.linspace(0.0, 2.0, N)k = t[1] - t[0]def odefunc(u, t):    ux = np.zeros(x.shape)    u[0] = 1 # boundary condition    for i in range(1,N-1):        ux[i] = (u[i+1] - u[i-1])/(2*h)         ux[i] -= k*(u[i+1] - 2*u[i] + u[i-1])/(2*h**2)    dudt = -ux    return dudt(note that I've renamed tspan to t for simplicity), you should get much better results.