Take a look at the simple fractional differential equation (FDE Initial Value Problem): $$_0D_x^{\frac{1}{2}}y(x)+2y(x)=0,~~x>0,~\text{and}~_0D_x^{-\frac{1}{2}}y(0)=1.$$ I solved the IVP using Laplace transform to get the exact solution:$$y(x)=x^{-\frac{1}{2}}E_{\frac{1}{2},\frac{1}{2}}(-2x^\frac{1}{2}),$$ where, ${\displaystyle E_{\alpha ,\beta }(z)=\sum _{k=0}^{\infty }{\frac {z^{k}}{\Gamma (\alpha k+\beta )}},}$ the Mittag-Leffler function.
As a student of numerical FDE, I am incredibly interested in comparing my numerical solution, and I tried to approximate $$_0D_x^{\frac{1}{2}}y(x)=\frac{d}{dx}\frac{1}{\sqrt \pi} \int_0^x(x-t)^{-\frac{1}{2}}y(t)dt=\frac{1}{\sqrt \pi}\frac{d}{dx}U(x),$$ by the central difference to follow, $$\frac{1}{\sqrt \pi}\frac{U(x_{i+1})-U(x_{i-1})}{2 \Delta x}+2y(x_i)=0,~\Delta x=x_{i+1}-x_i,$$ $$\Rightarrow \frac{1}{\sqrt \pi} \frac{\int_0^{x_{i+1}}(x_{i+1}-t)^{-\frac{1}{2}}y(t)dt-\int_0^{x_{i}}(x_{i}-t)^{-\frac{1}{2}}y(t)dt}{2\Delta x}+2y(x_i)=0,$$ $$\Rightarrow \int_0^{x_{i+1}}(x_{i+1}-t)^{-\frac{1}{2}}y(t)dt-\int_0^{x_{i}}(x_{i}-t)^{-\frac{1}{2}}y(t)dt+4\sqrt \pi\Delta x y(x_i)=0,$$ In the first look, I excited to club these integrals but it was a disaster due to the difference in the indices. Moreover, I am pretty scared about incorporating the initial condition. How can we proceed this central difference approximation.
Thanks in advance for any sort of help in the numerical analysis.