In this paper, we construct a complex reproducing kernel space for singular multi-point BVPs, and skillfully obtain reproducing kernel expressions. Then, we transform the problem into an equivalent operator equation, and give a numerical algorithm to provide the approximate solution. The uniform convergence of this algorithm is proved, and complexity analysis is done. Lastly, we show the validity and feasibility of the numerical algorithm by two numerical examples. Differential equations arise from various practical problems in mathematics and physics such as gas dynamics, nuclear physics, chemical reaction and geological prospecting etc. Multi-point BVPs can solve the contradictions in the process of actual research status efficiently and enhance their compatibility. Therefore, this problem has received a lot of attention of researchers. In [1 - 4] , the existence and uniqueness are investigated. In [5] , the author give the approximate solutions of a certain class of singular two-point (BVPs) by the Sinc-Galerkin method and homotopy-perturbation method. In [6] , the author solve two-point BVPs by variational iteration method. In [7] , the author present a method for solving a class of singular two-point BVPs based on cubic splines. For more information, please refer to [8 - 13] . Recently, in [14] , the author adopt differential transformation method to solve the two-point BVPs This method is able to provide an approximate solution of (1.1), while it has several disadvantages. This method is based on Taylor series which requires a high degree of smoothness. It has a local convergence region and the function f(x) are all polynomial functions. Now, we present a new algorithm to make up the deficiencies in [14] . It can be extended to singular multi-point BVPs. We take the problem where p(x), q(x), f(x) are continuous complex functions on (0, 1] and x = 0 is the singular value point of p(x), c, δ are complex constants and c ∈ (0, 1), δ ≠ 1. In this paper, we construct a complex reproducing kernel space [0, 1]. Then (1.2) can be transformed into an equivalent operator equation. Its approximate solution is provided. We also analyze the convergence and complexity of this algorithm. Finally, we give some numerical examples to verify the e®ectiveness of our algorithm. 2 .1. The complex reproducing kernel space [0, 1]. We define the inner product space [0, 1] = { u(x) } u ″ is absolutely continuous complex function, u ⁽³⁾ ∈ L ² [0, 1], u ′(0) = 0, u (1) = δu (c)}. Lemma 2.1. The space [0, 1] is a complex reproducing kernel space . The proof can be found in [15] . Next, we give the reproducing kernel function R_y(x) of [0, 1]. For each y ∈ [0, 1] and each u(x) ∈ [0, 1], by applying (2.1), we have that is where We can obtain Now, we have Since R_y(x) ∈ C ² , we get Similarly, if 0 ≤ y ≤ c , the function R_y(x) should satisfy the differential equation: and if c ≤ y ≤ 1, the function R_y(x) should satisfy the differential equation: If x ≠ y , (2.5) and (2.6) become Its characteristic equation is λ⁶ = 0, the characteristic roots are λ _i = 0, ( i = 1, 2, ⋯, 6). So, we assume that Since we get and In addition, we also need Combining(2.3), (2.4), (2.8) – (2.10) as well as y ∈ (0, c ) and y ∈ ( c , 1), we get 36 equations. We can find the undetermined coeffcients c_ij of (2.7) by solving the equations. If c = 1=2 and δ = 1 + 2 i , R_y(x) is the following expression 2.2. The complex reproducing kernel space [0, 1]. [0, 1] = { u(x) | u is absolutely continuous complex function, u ′ ∈ L ² [0, 1]}. The inner product is given by It is easy to prove that [0, 1] is a complex reproducing kernel space and its reproducing kernel is In this section, we investigate how to obtain approximate solutions of (1.2). First, we transform (1.2) into an equivalent operator equation (3.1). Then we give its approximate solution. Also, the convergence and complexity analysis are provided. 3.1. Equivalent operator equation. The equation (1.2) can be transformed into the following form: where α is constant and satisfy (b ≠ 0). Define linear operator : [0,1] → [0, 1] by Obviously, operator is bounded. The equation (1.2) can be converted into an equivalent operator equation: where f₁(x) = x^af(x) . 3.2. The numerical solution for operator equation (3.1). We choose a countable dense subset ⊂ (0, 1] and define 𝜓 _i(x) as Theorem 3.1. The function system is a complete system in the space [0, 1]. Proof . For an arbitrary i , we have, 0 = ⟨u(x), 𝜓_i(x)⟩ = ⟨u(x), (R_x(·))(x_i)⟩ = (⟨u(x),R_x(·)⟩)(x_i) = (u(·))(x_i) = (u)(x_i). Note that is dense in [0, 1], so ( u )( x ) = 0. By the existence of ⁻¹ , it follows that u ≡ 0. Therefore, is a complete system in [0, 1]. □ Furthermore, we obtain an orthogonal system of [0, 1] derived from Gram-Schmidt orthonormalization process from : Theorem 3.2. If is dense on [0, 1], then the solution of (1.2) is Proof . We expand u(x) into a Fourier series as follows Now, we can get the approximate solution u_n(x) by truncating the nth – term of the exact solution u(x) , 3.3. Theoretical analysis for our algorithm. Theorem 3.3. An approximate solution un(x) is uniform convergence to u(x) on [0, 1]. Moreover , , are both uniform convergence to u′(x) and u″(x) on [0, 1] . Proof . Note that u_n(x) = ⟨u_n, R_x⟩, u(x) = ⟨u,R_x⟩, and By applying Schwarz ’ s inequality and the boundedness of ( i = 0, 1, 2), we have So Theorem 3.4. The time complexity of the algorithm is O ( n ³ ). Proof . There are three steps to calculate the approximate solution u_n(x) of (1.2) . (1) Assume the number of multiplications required is C in one calculation of the inner product ⟨ φ_i,φ_i ⟩, then the total number of multiplications required is n ( n + 1) C =2 in calculation of all inner products. (2) Orthogonalization of the system needs 3 layers of nested loops, that is, the number of multiplication is (3) The number of multiplication is n ² when calculating u_n(x) using (3.2). To sum up, the total number of multiplication is In this section, some numerical examples are studied to demonstrate the accuracy of the present algorithm. Results obtained by this algorithm are compared with the exact solution of each example and are shown to be in good agreement with the exact solution. Example 1. Consider equation where f(x) = [(8+4 i )−(16−8 i ) −8 ie −(1−5 i ) e^x +(3−2 i ) x −((2−4i) + (4 + 8 i ) − 4 e ) x ² + . Its exact solution is Applying our algorithm and taking the number of nodes as n=50 and 100, the absolute errors of real part (a.e.Re) and the absolute errors of imaginary part (a.e.Im) are shown in Table 1 . It shows that the approximate solution is getting more and more accurate as n increases. Example 2. Consider equation where The function f(x) is continuous at x = , but not differentiable. So the method of [14] is invalid for example 2. While using our algorithm, we choose 100 points in (0, 1]. The numerical results | u ₁₀₀ − f | are given in the following Table 2 . In this paper, we present a new numerical algorithm in complex reproducing kernel space for singular multi-point BVPs. We give the rigorous theoretical analysis, the uniform convergence of the approximate solution. The numerical examples show that by using this algorithm we obtain better solution and fix the deficiencies of [14] .