## 3 Derivation of The Method  
### Definition of the Piecewise Spline Functions  
In this subsection, we use the theory of the analytic solution to construct the spline functions. Firstly, the interval \([a,b]\) is equally divided into \(n\) subintervals \([x_{i-1},x_{i}](i=1,2,...,n)\) such that \(\triangle:a=x_{0}<x_{1}<\cdots<x_{n}=b\) with \(h=x_{i}-x_{i-1}=(b-a)/n\), and we denote \(x_{i\pm 1/2}=x_{i}\pm\frac{1}{2}h\). If \(h\) is small, the value of \(x\) varies little at the interval \([x_{i-1},x_{i}]\), then \(P(x)\) and \(Q(x)\) can be seem as constants at the interval \([x_{i-1},x_{i}]\), therefore Equation (1) can also be seem as a second order linear differential equation with constant coefficients at the interval \([x_{i-1},x_{i}]\). We assume that the value of \(x\) at the interval \([x_{i-1},x_{i}]\) identically equals to the midpoint \(x_{i-1/2}\), hence Equation (1) can be seem as  
\[y^{\prime\prime}+P_{i-1/2}y=Q_{i-1/2}, \tag{4}\]  
where \(P_{i-1/2}=P(x_{i-1/2})\), \(Q_{i-1/2}=Q(x_{i-1/2})\). The roots of the characteristic equation of (4) is \(\lambda_{i}=\pm\sqrt{-P_{i-1/2}}\). Therefore, the solution of (4) can be obtained in the following three cases. **case 1. ** If \(P_{i-1/2}<0\), the characteristic roots \(\lambda_{i}\) are nonzero real numbers, it follows from Lemma 1 that \(y_{h}=a_{i}e^{\lambda_{i}(x-x_{i})}+b_{i}e^{-\lambda_{i}(x-x_{i})}\). Additionally, since the right side of Equation (4) identically equals to the constant \(Q_{i-1/2}\), we can assume that \(y_{p}=c\) according to Lemma 2. By substituting \(y_{p}\) into Equation (4), we have \(c=\frac{Q_{i-1/2}}{P_{i-1/2}}\). Therefore, the approximate function at the interval \([x_{i-1},x_{i}]\) can be defined as  
\[\tilde{y}_{i}(x)=a_{i}e^{\lambda_{i}(x-x_{i-1})}+b_{i}e^{-\lambda_{i}(x-x_{i-1 })}+\frac{Q_{i-1/2}}{P_{i-1/2}}.\]  
**case 2. ** If \(P_{i-1/2}>0\), then the characteristic roots \(\lambda_{i}\) are imaginary numbers, it follows from Lemma 1 that \(y_{h}=c_{i}\sin\lambda_{i}(x-x_{i-1})+d_{i}\cos\lambda_{i}(x-x_{i-1})\). Similarly to that described in Case 1, the particular solution is \(y_{p}=\frac{Q_{i-1/2}}{P_{i-1/2}}\). Therefore, the approximate function at the interval \([x_{i},x_{i+1}]\) can be defined as  
\[\tilde{y}_{i}(x)=c_{i}\sin\lambda_{i}(x-x_{i-1})+d_{i}\cos\lambda_{i}(x-x_{i- 1})+\frac{Q_{i-1/2}}{P_{i-1/2}}.\]  
**case 3. ** If \(P_{i-1/2}=0\), then \(\lambda_{i}=0\), the general solution of the homogeneous is \(y_{h}=e_{i}+f_{i}(x-x_{i-1})\). There are many particular solutions, here we take the simplest form \(y_{p}=g_{i}(x-x_{i-1})^{2}\). By substituting \(y_{p}\) into Equation (4), we yield \(g_{i}=\frac{1}{2}Q_{i-1/2}\). Therefore, the approximate function at the interval \([x_{i},x_{i+1}]\) can be defined as  
\[\tilde{y}_{i}(x)=\frac{1}{2}Q_{i-1/2}(x-x_{i-1})^{2}+e_{i}(x-x_{i-1})+f_{i},\]  
To summarize, the approximate function at the interval \([x_{i-1},x_{i}]\) is defined as  
\[\tilde{y}_{i}(x)=\left\{ \begin{array}{ll}a_{i}e^{\lambda_{i}(x-x_{i-1})}+b_ {i}e^{-\lambda_{i}(x-x_{i-1})}+\frac{Q_{i-1/2}}{P_{i-1/2}},&\text{if $P_{i-1/2}<0$;}\\ c_{i}\sin\lambda_{i}(x-x_{i-1})+d_{i}\cos\lambda_{i}(x-x_{i-1})+\frac{Q_{i-1/2} }{P_{i-1/2}},&\text{if $P_{i-1/2}>0$;}\\ \frac{1}{2}Q_{i-1/2}(x-x_{i-1})^{2}+e_{i}(x-x_{i-1})+f_{i},&\text{if $P_{i-1/2}=0$,} \end{array} \right. \tag{5}\]  
where \(a_{i},b_{i},c_{i},d_{i},e_{i}\) and \(f_{i}\) are the coefficients to be determined. Let the spline Function (5) satisfy the interpolation conditions  
\[\tilde{y}_{i}(x_{i-1})=y_{i-1},\tilde{y}_{i}(x_{i})=y_{i}. \tag{6}\]  
6.