[Python] Simple 1 D dispersion equation numerical solution

I am new to coding and trying to solve a simple 1D dispersion equation.

The equation and boundary conditions:

adC/dx = bd^2C/dx^2

x = hj, C = C0; x = - inf, C =0

The analytical solution is

C = C0 * exp (a/b*(x-hj))

Here is the code I have :

import numpy as np
import matplotlib.pyplot as plt

C0 = 3.5 # Concentration at injection point
h_j = 0.7 # Injection point
L = -100 # Location of closed boundary

alpha = 2.597
beta = 1.7

N=10000 # number of segments
x=np.linspace(L,h_j,N+1) # grid points
# determine dx
dx = x[1]-x[0]

# diagonal elements
a1 = -(alpha*dx)/(2*beta)
a2 = 2
a3 = alpha*dx/(2*beta)-1

# construct A
A = np.zeros((N+2, N+2))
A[0,0] = 1
A[N+1,N+1] = 1
for i in np.arange(1,N+1):
A[i,i-1:i+2] = [a1,a2,a3]
# construct B
B = np.zeros((N+2))
B[0] = 10^(-10)
B[N+1] = C0

xs=np.linspace(0,h_j,100)
Exact = C0*np.exp(alpha/beta*(xs-h_j))

#solve linear system Au=B

u = np.linalg.solve(A,B)

fig = plt.figure()
ax = fig.add_subplot(111)
res = ax.imshow(A, cmap = "jet", interpolation = 'nearest')
cb = fig.colorbar(res)

# drop the frist point at x_0, as it is outside the domain
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(x, u[1:N+2], label= "numerical solution")
ax.plot(xs,Exact, label= 'Analytical solution')
ax.legend()
ax.set_xlim(0,h_j)
ax.set_xlabel("z")
ax.set_ylabel("C(z)")
plt.show()




But the solved numerical solution is not the same as the analytical one. I cannot find the mistake in the codes or maybe in my math.

Continue reading...