I think I have found out the method to solve my question by the following codes

```
expr = fem.Expression(grad(u1), V1.element.interpolation_points())
w = Function(V)
u1_dx = w.sub(0).collapse()
u1_dx.interpolate(expr)
```

The whole example codes are be changed as follows:

```
'''
coordinate:
x->[0,10]
equations:
-u1'' - 4* u2'*u1' + 2*a*u1- 8*u1**3 + 6*u1**5 =0 , where the ' = d/dx,
# -u2''- 2/3 *u1'**2 =0 splited into
u3' + 2/3 *u1'**2 =0
u3 = u2'
boundary conditions
u1(0) =0, # u1'(10) =0
u2(0) =0, #u2'(0) =0
u3(0) =0
'''
import numpy as np
import ufl
from dolfinx import mesh, fem,log
from dolfinx.fem import Function, FunctionSpace
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from ufl import dx, split, TestFunctions,grad
from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
# parmeters
r_0 = 0.0
r_b = 10.0 # the coordinate of right boundary
a=0.2
b=2.0
c=1.0
beta=0.2
num =1000
msh = mesh.create_interval(comm=MPI.COMM_WORLD, points=(r_0, r_b), nx=num-1)
fdim = msh.topology.dim - 1
P1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)
V = FunctionSpace(msh, ufl.MixedElement(P1, P1,P1))
# the left boundary
facets_l = mesh.locate_entities_boundary(msh, dim=fdim, marker=lambda x: np.isclose(x[0], r_0))
dofs_u1 = fem.locate_dofs_topological(V=V.sub(0), entity_dim=fdim, entities=facets_l)
dofs_u2 = fem.locate_dofs_topological(V=V.sub(1), entity_dim=fdim, entities=facets_l)
dofs_u3 = fem.locate_dofs_topological(V=V.sub(2), entity_dim=fdim, entities=facets_l)
# the Dirichleb boundary conditions
bc_u1 = fem.dirichletbc(value=ScalarType(0), dofs=dofs_u1, V=V.sub(0)) # u1 = phi
bc_u2 = fem.dirichletbc(value=ScalarType(0), dofs=dofs_u2, V=V.sub(1)) # u2 = A
bc_u3 = fem.dirichletbc(value=ScalarType(0), dofs=dofs_u3, V=V.sub(2)) # u2 = A
bcs = [bc_u1, bc_u2, bc_u3]
# Define functions
u = Function(V)
# Split system functions to access components
#phi, m, sigma = u1, u2, u3
u1, u2,u3 = split(u)
# Define test functions
v1, v2, v3 = TestFunctions(V)
# Define variational problem
F1 = v1.dx(0)*u1.dx(0) + (-4.0 * u2.dx(0)*u1.dx(0) +2.0*a*u1- 4.0*b*u1**3 + 6.0*c*u1**5)*v1
F2 = (u3.dx(0)+ 2.0/3.0*u1.dx(0)**2)*v2
F3 = (u3-u2.dx(0))*v3
F = (F1+F2+F3)*dx
# guess solutions
x=np.linspace(r_0,r_b,num)
y1=[0,1,1,1,1,1]
y3=[0,0,-0.8,-1.2,-1.6,-2.0]
x1=np.linspace(r_0,r_b,num=6)
f1=interp1d(x1,y1,kind='linear')
f3=interp1d(x1,y3,kind='linear')
y1_ = f1(x)
y2_ = np.diff(y1_)
y2_ = np.insert(y2_,0,0)
y3_ = f3(x)
#yguess=guess(x)
V1, dofs0 = V.sub(0).collapse()
V2, dofs1 = V.sub(1).collapse()
V3, dofs2 = V.sub(2).collapse()
u.x.array[dofs0]=y1_ # yguess[0]
u.x.array[dofs1]=y2_ # yguess[1]
u.x.array[dofs2]=y3_ # yguess[2]
# Create nonlinear problem and Newton solver
problem = NonlinearProblem(F, u, bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-8
solver.nonzero_initial_guess = True
solver.max_it = 100
solver.report = True
out = solver.solve(u)
u1_ = u.x.array[dofs0]
u2_ = u.x.array[dofs1]
u3_ = u.x.array[dofs2]
# u1_ = u.sub(0).collapse()
expr = fem.Expression(grad(u1), V1.element.interpolation_points())
w = Function(V)
u1_dx = w.sub(0).collapse()
u1_dx.interpolate(expr)
expr2 = fem.Expression(grad(u2), V2.element.interpolation_points())
u2_dx = w.sub(1).collapse()
u2_dx.interpolate(expr2)
plt.figure(1)
plt.plot (x, u1_,label="u1")
plt.plot (x, u2_,label="u2")
plt.plot (x, u3_,label="u3") # the solution of u3 is the u2_dx theoritical
plt.plot (x, (u3_ - u2_dx.x.array[:]), label="u3-u2_dx")
plt.plot (x, u1_dx.x.array[:],label="u1_dx")
plt.plot (x, u2_dx.x.array[:],label="u2_dx")
plt.legend()
plt.show()
plt.close ( )
```

The output fig is