I want to solve the equation: $$ u_t +\partial_x\left(e^{-x}u \right) = 0$$ with the BC $u(0,z,t)=u_0$ on the domain $D=[0,10]\times[0,1]$. My attempt is to use the following approximations: $$ \frac{u^{n+1}_{j,k}-u^n_{j,k}}{\Delta t} + e^{-x_j} \frac{u^{n+1}_{j,k}-u^{n+1}_{j-1,k}}{\Delta x} - e^{-x_j} u^{n+1}_{j,k}=0 $$ However, I wonder what the order of the scheme in $\Delta x$ will be. I think the upwind scheme should be more accurate than the centered difference scheme since it's an advection-dominated problem but my code doesn't seem to be working correctly and I'm not sure why. Any ideas on how I can increase the accuracy or should this be sufficient?
EDIT
My code results in the following solution at $t_F=10$:
Doubling the number of points in $x$ for each simulation and comparing to the solution from the most refined grid yields a convergence plot as shown below:
Even further, I can use the solution data to approximate terms in the PDE and so I can actually obtain a value for the LHS of the PDE at each point in my discretized grid. Doing this yields:
So the convergence is exactly as expected (actually slightly better), the error as measured by the solution is quite low except near the transition region (where it is expected to be highest), yet the numerical solution looks nothing like the exact, analytical solution. So my question is why? Is there an error in my code?
For reference, here I used the initial condition: $$ u(x,0) = u_0 - \frac{u_0}{1+e^{-4(x-1)}} $$ where $u_0\in [0,1]$ is a constant.