Time dependent problem, Converges to one point for some reason

General
1

Why does the graph converge at different points in time to the same point? When “xi” is responsible for the spatial coordinate. The graph should with the right border tends to 0.

from dolfin import *
import matplotlib.pyplot as plt
import numpy as np
import scipy.special
import math
import os

plt.figure(figsize=(10, 5), dpi=80)

set_log_level(LogLevel.ERROR)
alpha = 1

T = 1
M = 100
dt = T/M
nx = 65
t = -1

k0 = 1
A = 1
sigma = 2
t_f = 0
L_f = np.sqrt(2*k0*(A**sigma)*(sigma+2)/sigma)
print('L_f =', L_f)

n = -1./(0.5*sigma) 

xi = Expression('x[0]/(pow(k0, 1/2)*pow(A, sigma/2)*pow(-t, (1+n*sigma)/2))', degree=10, sigma=sigma, n=n, t=t, k0=k0, A=A)

u0 = Expression('A*pow(-t,n)*pow((1.-xi/(pow(2*k0*(pow(A,sigma))*(sigma+2)/sigma, 0.5))),(2./sigma))', degree=10, sigma=sigma, t=t, xi=xi, n=n, A=A, k0=k0) #точное решение

mesh = IntervalMesh(nx, 0, L_f)
V = FunctionSpace(mesh, 'CG', 2)

def b1(x, on_boundary):
    if x[0] <= L_f:
    	return near(x[0], 0)
    if x[0] > L_f:
    	return 0
def b2(x, on_boundary):
    if x[0] <= L_f:
    	return near(x[0], L_f)
    if x[0] > L_f:
    	return 0
bc1 = DirichletBC(V, u0, b1)
bc2 = DirichletBC(V, 0, b2)

v = TestFunction(V)
u = Function(V)

def q(uk):
	return uk**sigma

ut = [interpolate(u0, V) for j in range(M)]

for k in range(1,M):
    t = t + dt   
    xi_values = xi.compute_vertex_values(mesh)
    xis = xi(1,66)
    xi.t = t
    u0.t = t
    u0.xi = xi
    
    F = u*v*dx + (dt**alpha)*dot(q(ut[k-1])*grad(u), grad(v))*dx - (ut[0] - sum( (-1)**j * scipy.special.binom(alpha, j) * (ut[k-j]-ut[0]) for j in range(1, k)) )*v*dx
    
    J = derivative(F, u)
    problem = NonlinearVariationalProblem(F, u, [bc1,bc2], J)
    solver  = NonlinearVariationalSolver(problem)
    prm = solver.parameters
    prm['newton_solver']['absolute_tolerance'] = 1E-5
    prm['newton_solver']['relative_tolerance'] = 1E-4
    prm['newton_solver']['maximum_iterations'] = 50000
    prm['newton_solver']['relaxation_parameter'] = 1.0
    solver.solve()

    ut[k].assign(u)

    error_L2 = errornorm(u0, u, 'L2')
    error_H1 = errornorm(u0, u, 'H1')

    plot(u, label='t=%.2f'%(t))

print(' ')
print (f'alpha={alpha}, \u0394t={dt}, \u0394x={1/nx}, M={M}, error_H1={error_H1}, error_L2={error_L2}')

plt.title(f'alpha={alpha}, t_f={t_f}, \u0394x={1/nx}')

plt.grid(True)
plt.ylim(-0.5, 25)
plt.show()

print('complete')

p