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.

You have a number of errors.

Your matrix coefficients are wrong - try discretising those first and second derivatives again carefully.

Your boundary values are fixed (Dirichlet BC). So you only have N-1 varying nodes. So your matrix should be (N-1)x(N-1) not (N+1)x(N+1). The boundary conditions are effectively transferred to be the only non-zero source terms (RHS).

Your leftmost boundary condition might as well be 0; then you don't have to worry about how you incorrectly wrote 1e-10. Everything is forced from the right boundary - as it would be physically.

Finally, you are using "central differencing" for the advection term (the first derivative on the LHS). This is dubious physically (information goes in the direction of flow) and may, for coarse grids, lead to unphysical numerical wiggles. In the long run, consider "upwind"/one-sided differencing. In the present case, your N is large enough, and dx small enough, not to cause a problem. However, if you, say, multiply alpha by 1000 you will get an ... interesting ... result.

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=5000 # 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) - 1
a2 = 2
a3 =  alpha*dx/(2*beta) - 1

# construct A
A = np.zeros( (N-1,N-1) )    # discretised derivative matrix
B = np.zeros( (N-1) )        # RHS (source term)
u = np.zeros( (N+1) )        # solution vector (NOTE: includes boundary nodes)
# Interior grid nodes
for i in range( 1, N-2 ):
    A[i,i-1] = a1
    A[i,i  ] = a2
    A[i,i+1] = a3
# Boundary nodes
A[0  ,0  ] = a2;   A[0  ,1  ] = a3;    B[0  ] = -a1 * 0 ;    u[0] = 0       # left  boundary
A[N-2,N-3] = a1;   A[N-2,N-2] = a2;    B[N-2] = -a3 * C0;    u[N] = C0      # right boundary
# Solve linear system Au=B (NOTE: interior nodes only; boundary values are Dirichlet BC)
u[1:N] = np.linalg.solve(A,B)

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

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

fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(x, u    , 'bo', label= 'numerical solution'  )
ax.plot(xs,Exact, 'r-', label= 'Analytical solution' )
ax.legend()
ax.set_xlim(0,h_j)
ax.set_xlabel("z")
ax.set_ylabel("C(z)")
plt.show()