Numerical solution for a non-linear Fractional Differential Equation (FDE)

Question

As shown below, a neat explicit expression is obtained for F=2, however an exact solution is not present for 1< F < 2.

How do we obtain numerical values for F = 1.5 (for instance)?

There have been suggestions of approximate methods (Adomian decomposition method (ADM), the variational iteration method (VIM)), these could involve time-consuming data crunching, and not sure if such techniques are already utilized in Mathematica.

The exact solution for the F=2 case does allow the testing of these possible alternative methods, hopefully one where rapid convergence is noted. An algorithm based on the Adomian decomposition method and the Padé approximants technique for nonlinear fractional differential equations with initial or boundary conditions is described here, but I am unable to access this code.

F = 2;


sol = DSolve[{CaputoD[y[x], {x, F}] - 9/4 Sqrt[y[x]] - y[x] == 0, 
   y[0] == 1, y'[0] == 2}, y[x], x]


Table[{x, Evaluate[y[x] /. sol[[1]]]}, {x, 0, 1, 0.1}]

Plot[Evaluate[y[x]  . sol[[1]]], {x, 0, 1}]

> (* {{y[x] ->     1/16 E^(-2 x) (E^x + 54 E^(2 x) + 81 E^(3 x) - 
>       12 Sqrt[E^(3 x) (1 + 9 E^x)^2])}} *)
> 
> (* {{0., 1.}, {0.1, 1.21698}, {0.2, 1.47099}, {0.3, 1.76705}, {0.4, 
>   2.11074}, {0.5, 2.50829}, {0.6, 2.96662}, {0.7, 3.49344}, {0.8, 
>   4.09733}, {0.9, 4.78782}, {1., 5.57553}} *)

Answer 1

There are several numerical methods to solve fractional DE and PDE, including predictor-corrector method described here, and wavelets method discussed in our paper and here. First we solve this problem for $1\le F\le 2$ with Haar wavelets. Note, that in this case we have the second order problem, and it is why we need two boundary conditions. In a case $F=2$ we have CaputoD[y[x], {x, F}]=D[y[x], {x, 2}], therefore problem turns to second order ODE with known solution. This solution we gonna use as a test. With Haar wavelets we can compute CaputoD[y[x], {x, q}] as follows

h[x_, k_, m_] := 
  WaveletPsi[HaarWavelet[], m x - k, WorkingPrecision -> Infinity];

p[x_, k_, m_] := 
 Piecewise[{{(1 + k - m*x)/m, k >= 0 && 1/m + (2*k)/m - 2*x < 0 && 
           1/m + k/m - x >= 0 && m > 0}, {(-k + m*x)/m, 
    k >= 0 && 1/m + (2*k)/m - 2*x >= 0 && 
           k/m - x < 0 && 1/m + k/m - x >= 0 && m > 0}}, 0];
h1[x_] := WaveletPhi[HaarWavelet[], x, WorkingPrecision -> Infinity];

p1[x_] := Piecewise[{{1, x > 1}}, x];

pc[t_, k_, m_, q_] := 
 Piecewise[{{-(t^(1 - q)/(-1 + q)), k == 0 && 1/m - 2*t >= 0 && 
          m > 0 && t > 0 && 1/m - t >= 0}, 
      {-((m^(-1 + q)*(1/(-k + m*t))^(-1 + q))/(-1 + q)), 
        k > 0 && 1/m + (2*k)/m - 2*t > 0 && k/m - t < 0 && m > 0 && 
          1/m + k/m - t > 0}, 
      {(-t^q + 2*m*t^(1 + q) - m*t*(-(1/(2*m)) + t)^q)/
          (t^q*(-(1/(2*m)) + t)^q*(m*(-1 + q))), 
        k == 0 && m > 0 && 1/m - 2*t < 0 && 1/m - t >= 0}, 
      {(1/(-1 + q))*((2^(-1 + q)*m^(-1 + 2*q)*(-(-(k/m) + t)^q - 
                   2*k*(-(k/m) + t)^q + 2*m*t*(-(k/m) + t)^q + 
                   2*k*(-((1/2 + k)/m) + t)^q - 
           2*m*t*(-((1/2 + k)/m) + t)^
                       q))/((1 + 2*k - 2*m*t)*(k - m*t))^q), 
        k > 0 && 1/m + (2*k)/m - 2*t == 0 && m > 0 && 
          1/m + k/m - t > 0}, 
      {-((1/(-1 + q))*((2^(-1 + q)*m^(-1 + 2*q)*
                   (-2*(-((1/2 + k)/m) + t)^
               q*((1 + 2*k - 2*m*t)*(k - m*t))^
                          q - 2*k*(-((1/2 + k)/m) + t)^q*
                        ((1 + 2*k - 2*m*t)*(k - m*t))^q + 
                      
             2*m*t*(-((1/2 + k)/m) + t)^q*((1 + 2*k - 2*m*t)*
                             (k - m*t))^q + (-((1 + k)/m) + t)^q*
                        ((1 + 2*k - 2*m*t)*(k - m*t))^q + 
                      
             2*k*(-((1 + k)/m) + t)^q*((1 + 2*k - 2*m*t)*(k - m*t))^
                          q - 2*m*t*(-((1 + k)/m) + t)^q*
                        ((1 + 2*k - 2*m*t)*(k - m*t))^
               q + (-(k/m) + t)^q*
                        ((1 + 2*k - 2*m*t)*(1 + k - m*t))^q + 
                      
             2*k*(-(k/m) + t)^q*((1 + 2*k - 2*m*t)*(1 + k - m*t))^
               q - 
                      
             2*m*t*(-(k/m) + t)^q*((1 + 2*k - 2*m*t)*(1 + k - m*t))^
                          q - 2*k*(-((1/2 + k)/m) + t)^q*
                        ((1 + 2*k - 2*m*t)*(1 + k - m*t))^q + 
                      
             2*m*t*(-((1/2 + k)/m) + t)^q*((1 + 2*k - 2*m*t)*
                             (1 + k - m*t))^
               q))/(((1 + 2*k - 2*m*t)*(k - m*t))^q*
                   ((1 + 2*k - 2*m*t)*(1 + k - m*t))^q))), 
        k > 0 && m > 0 && 1/m + (2*k)/m - 2*t <= 0 && 
          1/m + k/m - t <= 0}, 
      {-((1/(2*m*(-1 + q)))*((2^q*m^(2*q)*t^q*(-(1/m) + t)^q*
                     (-(1/(2*m)) + t)^q - 
           2^(1 + q)*m^(1 + 2*q)*t^(1 + q)*
                     (-(1/m) + t)^q*(-(1/(2*m)) + t)^q - 
           2^(1 + q)*m^(2*q)*
                     t^q*(-(1/(2*m)) + t)^(2*q) + 
           2^(1 + q)*m^(1 + 2*q)*
                     t^(1 + q)*(-(1/(2*m)) + t)^(2*q) + 
                   t^q*((-1 + m*t)*(-1 + 2*m*t))^q - 2*m*t^(1 + q)*
                     ((-1 + m*t)*(-1 + 2*m*t))^q + 
           2*m*t*(-(1/(2*m)) + t)^q*
                     ((-1 + m*t)*(-1 + 2*m*t))^q)/(t^
            q*(-(1/(2*m)) + t)^q*
                   ((-1 + m*t)*(-1 + 2*m*t))^q))), 
        k == 0 && 1/m - 2*t < 0 && 1/m - t < 0 && m > 0}, 
      {(1/(-1 + q))*((2^(-1 + q)*m^(-1 + q)*((-m^q)*(-(k/m) + t)^q - 
                   2*k*m^q*(-(k/m) + t)^q + 
           2*m^(1 + q)*t*(-(k/m) + t)^q + 
                   2*k*m^q*(-((1/2 + k)/m) + t)^q - 2*m^(1 + q)*t*
                     (-((1/2 + k)/m) + t)^
             q - ((1 + 2*k - 2*m*t)*(k - m*t))^q*
                     (1/(-1 - 2*k + 2*m*t))^q - 
                   2*k*((1 + 2*k - 2*m*t)*(k - m*t))^q*
                     (1/(-1 - 2*k + 2*m*t))^q + 
                   2*m*t*((1 + 2*k - 2*m*t)*(k - m*t))^q*
                     (1/(-1 - 2*k + 2*m*t))^q))/((1 + 2*k - 
            2*m*t)*(k - m*t))^
               q), 1/m + (2*k)/m - 2*t < 0 && k > 0 && m > 0 && 
          1/m + k/m - t > 0}}, 0]

pc1[t_, q_] := Piecewise[{{-(t^(1 - q)/(-1 + q)), t <= 1}}, 
    -(((-1 + t)^q*t + t^q - t^(1 + q))/((-1 + t)^q*t^q*(-1 + q)))];  

Here p, p1 are integrals of h,h1, and pc, pc1 are Caputo derivatives of h, h1. With these functions we can compute numerical solution for $0\le q\le 1$ by reducing order FDE as

eqs={CaputoD[z[x],{x,q}] - 9/4 Sqrt[y[x]] - y[x] == 0, z[x]==y'[x],
   y[0] == 1, z[0] == 2}

Finally we define numerical functions y[x],z[x] and first derivatives y1[x], z1[x], and compute eqs in collocation points xcol

J = 4; M = 2^J; dx = 1/(2*M); xl = Table[l dx, {l, 0, 2 M}]; xcol = 
 Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, 2 M + 1}]; 
z1[x_, q_] := 
 Sum[v[i, j] pc[x, i, 2^j, q], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
  v1 pc1[x, q]; 
z[x_] := 
 Sum[v[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
  v1 p1[x] + v0;
y1[x_] := 
 Sum[u[i, j] h[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
  u1 h1[x]; 
y[x_] := 
 Sum[u[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
  u1 p1[x] + u0;
varM = Join[{v0, v1, u0, u1}, 
  Flatten[Table[v[i, j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}]], 
  Flatten[Table[u[i, j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}]]]; 
eq[q_] := 
 Flatten[Table[{1/Gamma[1 - q] z1[x, q] - 9/4 Sqrt[Abs[y[x]]] - y[x] ==
      0, y1[x] - z[x] == 0}, {x, xcol}]]; bc = {y[0] == 1, z[0] == 2};

The system of algebraic equations can be solved with FindRoot for several q

qlist = {99/100, 1/2, 1/3};

Do[
  eqM[i] = Join[eq[i], bc];
  solv[i] = 
   FindRoot[eqM[i], Table[{varM[[j]], 1/10}, {j, Length[varM]}], 
    MaxIterations -> 1000, WorkingPrecision -> 30];, {i, 
   qlist}]  

Visualization together with test data

Do[lst[i] = 
   Table[{x , Evaluate[y[x] /. solv[i]]}, {x, 0, 1, .01}];, {i, qlist}]

test = {{0., 1.}, {0.1, 1.21698}, {0.2, 1.47099}, {0.3, 
    1.76705}, {0.4, 2.11074}, {0.5, 2.50829}, {0.6, 2.96662}, {0.7, 
    3.49344}, {0.8, 4.09733}, {0.9, 4.78782}, {1., 5.57553}};

Show[ListLinePlot[Table[lst[i], {i, qlist}], Frame -> True, 
  FrameLabel -> {"x", "y"}, 
     PlotRange -> All, 
  PlotLegends -> 
   Table[Row[{"F = ", qlist[[i]] + 1}], {i, Length[qlist]}]], 
 ListPlot[test, PlotStyle -> Red]] 

Predictor-corrector code is based on the pseudocode published in the end of paper

g[t_, y_] := 9/4 Sqrt[y] + y;

\[Alpha] = 199./100;
h = 1./100;
y[0] = y0[0] = 1;
y0[1] = 2; n = 100;
For[k = 1, k <= n, k++, b[k] = k^\[Alpha] - (k - 1)^\[Alpha]; 
  a[k] = -(2*k^(\[Alpha] + 1)) + (k - 1)^(\[Alpha] + 1) + (k + 
       1)^(\[Alpha] + 1);];
For[j = 1, j <= n, j++, 
  p[j] = (h^\[Alpha]*Sum[b[j - A]*g[A*h, y[A]], {A, 0, j - 1}])/
      Gamma[\[Alpha] + 1] + y0[0] + j h y0[1];
   
   y[j] = (h^\[Alpha]*(Sum[a[j - f]*g[f*h, y[f]], {f, 1, j - 1}] + 
          g[h*j, p[
            j]] + ((j - 1)^(\[Alpha] + 1) - (-\[Alpha] + j - 1)*
              j^\[Alpha])*g[0, y[0]]))/Gamma[\[Alpha] + 2] + y0[0] + 
     j h y0[1];
   ;];

Visualization together with test data

lst = Table[{j h, y[j]}, {j, n}];
test = {{0., 1.}, {0.1, 1.21698}, {0.2, 1.47099}, {0.3, 
   1.76705}, {0.4, 2.11074}, {0.5, 2.50829}, {0.6, 2.96662}, {0.7, 
   3.49344}, {0.8, 4.09733}, {0.9, 4.78782}, {1., 5.57553}}; Show[
 ListLinePlot[lst], ListPlot[test, PlotStyle -> Red]] 

Finally note, that solver proposed by xzczd is not so good compare to wavelets and predictor-corrector method and needs to be improved. Code is given by

F = 3/2; yfunc = 
 x |-> a + b x + 
   c x^2; ivp = {CaputoD[y[x], {x, F}] - 9/4 Sqrt[y[x]] - y[x] == 0 /. 
   x -> x0 + dx, y[x0] == y0, y'[x0] == dy0};
solrule = 
 Solve[ivp /. y -> yfunc, {a, b, 
   c}]; (Subtract @@@ ivp /. y -> yfunc /. # /. {x0 -> 0, y0 -> 1, 
      dy0 -> 2, dx -> 0.1, a -> 1, b -> 2} // N) & /@ solrule; cf = 
 Compile[{x0, y0, dy0, dx}, 
  Evaluate[{x0 + dx, yfunc[x0 + dx], yfunc'[x0 + dx], dx} /. 
    solrule[[-1]]]];

solver[x0_, y0_, dy0_, dx_, tend_] := 
 MapAt[List, 
    NestList[cf @@ # &, {x0, y0, dy0, dx}, 
       tend/dx // Round]\[Transpose] // Most, {1, All}] // Transpose //
   Interpolation
tend = 1;
nsol = solver[0, 1, 2, 0.01, tend]
nsol // ListPlot

In the picture below shown Figure 1 with wavelets method results, predictor-corrector method for F=3/2 - dashed red line and nsol - doted line. We see that two methods are in a good agreement while nsol is out of range.

Updated code by xzczd is in a good agreement with Haar wavelets method and predictor-corrector method as well. Code

n = 16;
range = Range[0, n - 1];
coef = c /@ range;
poly[x_] = coef . x^range;
CGLGrid[xl_, xr_, n_Integer /; n > 1] := 
 1/2 (xl + xr + (xl - xr) Cos[(\[Pi] Range[0, n - 1])/(n - 1)])
grid = CGLGrid[0, 1, n];

F = 3/2;
{eq, ic} = {CaputoD[y[x], {x, F}] - 9/4 Sqrt[y[x]] - y[x] == 
    0, {y[0] == 1, y'[0] == 2}};

approx = eq /. y -> poly; // AbsoluteTiming
cguess = 1;
crule = FindRoot[
    Table[approx, {x, grid[[3 ;;]]}]~
     Join~(ic /. y -> poly), {#, cguess} & /@ coef]; // AbsoluteTiming


pseudo = 
 Plot[poly[x] /. crule // Evaluate, {x, 0, 1}, PlotRange -> All, 
  GridLines -> Automatic, PlotStyle -> Dashed]

In the picture below the dashed line is the numerical solution computed with xzczd code compare to wavelets (left) and predictor-corrector method (right):

Answer 2

Update 2: A general-purpose IVP solver for FDE

The following is a general-purpose implementation for the predictor-corrector method in this paper linked by Alex:

BeginPackage["fdeivp`"];
fdesolve; fdesolvecore;
Begin["`private`"];
Clear[fdesolve, fdesolvecore];
Options[fdesolve] = {Compiled -> True};
fdesolve[α_, f_, y0_?NumericQ, rest__] := fdesolve[α, f, {y0}, rest]
fdesolve[α_, f_, y0_?VectorQ, rest__] /; Length@y0 != Ceiling@α := 
 fdesolve[α, f, {y0}, rest]
fdesolve[α_, f_, y0_, T_, n_Integer : 25, OptionsPattern[]] :=
 With[{size = Rest@Dimensions@y0},
  fdesolvecore[α, f, size, OptionValue@Compiled][y0, T, n]]
   
fdesolvecore[α_, f_, size_, tf : True | False] := 
  fdesolvecore[α, f, size, tf] = 
   With[{m = Ceiling[α], dim = 1 + Length@size}, 
    Hold@function[argu, with[{h = T/n,                 
               b = Table[If[k == 0, 0, k^α - (k - 1)^α], {k, 0, n}],                 
               a = Table[If[k == 0, 0, (k + 1)^(α + 1) - 
                             2 k^(α + 1) + (k - 1)^(α + 1)], {k, 0, n}]},
              module[{y = Table[0., {n + 1}, ##], p},
               y[0] = y0[0];                  
               Do[p = Sum[(j h)^k/k! y0[k], {k, 0, m - 1}] + 
                  h^α/Gamma[α + 1] Sum[b[j - k] f[k h, y[k]], {k, 0, j - 1}];   
                y[j] = Sum[(j h)^k/k! y0[k], {k, 0, m - 1}] + 
                  h^α/Gamma[α + 2] (f[j h, p] + ((j - 1)^(α + 1) - (j - 1 - α) j^α) f[
                    0, y[0]] + Sum[a[j - k] f[k h, y[k]], {k, j - 1}]),                
                {j, n}]; y]]] & @@ (List /@ size) /. (head : y | y0 | a | b)[i_] :> 
          head[[i + 1]] /. {module -> Module, with -> With} /. {function -> 
         If[tf, Compile, Function], 
        argu -> If[tf, {{y0, _Real, dim}, {T, _Real}, {n, _Integer}}, {y0, T, n}]} /. 
      Factorial[n_] :> Gamma[1 + n] // ReleaseHold];
End[];
EndPackage[];

With this package, OP's problem can be solved as follows:

order = 1.5; rhsfunc = {t, y} |-> 9/4 Sqrt[y] + y; iclst = {1, 2}; tend = 1;

sollst = fdesolve[order, rhsfunc, iclst, tend];

The output is a List, it can be visualized with e.g.

ListPlot[sollst, DataRange -> {0, tend}]

If you prefer a function relationship as the output, build a interpolating function with e.g.

solfunc = ListInterpolation[sollst, {0, tend}]

---

Update

Though my old answer is wrong, I manage to think out a correct one.

The idea is again simple: given that usually the solution of differential equation can be approximated (without Runge's phenomenon) with n-th order polynomial whose coefficient is determined by function values at Chebyshev-Gauss-Lobatto grid (notice this is essentially the spirit of pseudospectral method), and CaputoD is capable of calculating derivative of polynomial, we assume the solution of the initial value problem (IVP) in the whole interval $[0,1]$ can be approximated by

$$y=c_0+c_1 x+…+c_{n-1} x^{n-1}$$

n = 16;
range = Range[0, n - 1];
coef = c /@ range;
poly[x_] = coef . x^range;

Then substitute it into the differential equation and i.c.s, and calculate the specific algebraic equations at Chebyshev-Gauss-Lobatto grid:

CGLGrid[xl_, xr_, n_Integer /; n > 1] := 
 1/2 (xl + xr + (xl - xr) Cos[(π Range[0, n - 1])/(n - 1)])
grid = CGLGrid[0, 1, n];

F = 3/2;
{eq, ic} = {CaputoD[y[x], {x, F}] - 9/4 Sqrt[y[x]] - y[x] == 0, 
            {y[0] == 1, y'[0] == 2}};

{approxeq, approxic} = {eq, ic} /. y -> poly; // AbsoluteTiming
cguess = 1;
crule = FindRoot[
    Table[approxeq, {x, grid[[3 ;;]]}]~Join~approxic, {#, cguess} & /@ 
     coef]; // AbsoluteTiming

pseudo = Plot[poly[x] /. crule // Evaluate, {x, 0, 1}, 
     PlotRange -> All, GridLines -> Automatic]

Now the solution is consistent with those shown by Alex.

BTW, the following is a compiled version of the predictor-corrector method shown in Alex's answer:

g[t_, y_] := 9/4 Sqrt[y] + y;

alexsolver = Hold@Compile[{{n, _Integer}, α}, Module[{
         h = 1./n,
         y = Table[0., {n + 1}],
         y0 = Table[0., {n + 1}],
         a = Table[0., {n + 1}],
         b = Table[0., {n + 1}],
         p = Table[0., {n + 1}], k, j},
        y[0] = y0[0] = 1;
        y0[1] = 2;
        For[k = 1, k <= n, k++, b[k] = k^α - (k - 1)^α;
         a[
           k] = -(2*k^(α + 1)) + (k - 1)^(α + 1) + (k + 1)^(α + 
              1);];
        For[j = 1, j <= n, j++, 
         p[j] = (h^α*Sum[b[j - A]*g[A*h, y[A]], {A, 0, j - 1}])/
            Gamma[α + 1] + y0[0] + j h y0[1];
         y[
           j] = (h^α*(Sum[a[j - f]*g[f*h, y[f]], {f, 1, j - 1}] + 
                g[h*j, p[j]] + ((j - 1)^(α + 1) - (-α + j - 1)*
                    j^α)*g[0, y[0]]))/Gamma[α + 2] + y0[0] + j h y0[1]];
        y], CompilationTarget -> "C", RuntimeOptions -> "Speed"] /. 
     DownValues@g /. (head : y | y0 | a | b | p)[i_] :> head[[i + 1]] // ReleaseHold;

n = 10^4; h = 1/n; F = 3/2; ylst = alexsolver[n, F]; // AbsoluteTiming

alexsol = ListInterpolation[ylst, {0, 1}];
pseudo~Show~ListLinePlot[alexsol, PlotStyle -> {Red, Dashed}]

---

Old Incorrect Answer

The following answer is wrong, because definition of CaputoD involves integration from 0 to x, so we can not approximate CaputoD only with function value in a small interval $[x_0,x_0+dx]$ (except for $x_0=0$).

Since currently NDSolve cannot handle fractional differential equation, I'd like to add a primary but working home-made solution.

The idea is simple: given that usually the solution of differential equation can be approximated with piecewise polynomial (notice this is essentially the spirit of finite difference method (FDM)), and CaputoD is capable of calculating derivative of polynomial, we assume that the solution of the initial value problem (IVP) in a small interval $[x_0,x_0+dx]$ can be approximated by

$$y\approx a+bx+cx^2$$

yfunc = x |-> a + b x + c x^2;

Substitue this into the IVP and solve:

ivp = {CaputoD[y[x], {x, F}] - 9/4 Sqrt[y[x]] - y[x] == 0 /. x -> x0 + dx, 
   y[x0] == y0, y'[x0] == dy0};
solrule = Solve[ivp /. y -> yfunc, {a, b, c}]; // AbsoluteTiming 

Solve::nongen There may be values of the parameters for which some or all solutions are not valid.

(* {0.0378036, Null} *)

Hmm… something improper seems to happen in Solve. Further check via Dimensions@solrule suggests 2 solutions are found, which is a bit counterinituive. I don't bother to looking into what happens in Solve, so I simply check the solutions numerically:

(Subtract @@@ ivp /. y -> yfunc /. # /. {x0 -> 0, y0 -> 1, dy0 -> 2, 
      dx -> 0.1, a -> 1, b -> 2} // N) & /@ solrule
(* {{-4.89298, 0., 0.}, {-4.44089*10^-16, 0., 0.}} *)

OK, seems that only the last solution is valid. Extract it, and substitute it into yfunc, we obtain the approximate solution at $x_0+dx$. Needless to say, the process can be done repeatively, so we can obtain the numeric solution in arbitrary interval. The programming is straightforward:

cf = Compile[{x0, y0, dy0, dx}, 
  Evaluate[{x0 + dx, yfunc[x0 + dx], yfunc'[x0 + dx], dx} /. solrule[[-1]]]]

solver[x0_, y0_, dy0_, dx_, tend_] := 
 MapAt[List, 
    NestList[cf @@ # &, {x0, y0, dy0, dx}, tend/dx // Round]\[Transpose] // Most, 
       {1, All}] // Transpose // Interpolation

Here, 4th syntax in document of Interpolation is used to take the derivative into account when building InterpolatingFunction.

Usage of solver:

tend = 1;      
nsol = solver[0, 1, 2, 0.1, tend]
nsol // ListPlot

I've used an undocumented syntax of ListPlot here. See this post for more info.

Finally let's check the convergency of this method:

With[{lst = 10.^-Range[5]}, 
 Table[solver[0, 1, 2, dx, tend], {dx, lst}] // 
  ListLinePlot[#, PlotLegends -> lst] &]

Um… not too bad (at least in my view). As mentioned in the beginning, this is merely a primary solution. Seeing it working, I'm satisfied enough.

Answer 3

Good reference is provided by "A review of the Adomian decomposition method and its applications to fractional differential equations ", Commun.Frac.Calc.3 (2) (2012) 73-99

Note that the non-linear operator is taken to be semi-analytic, and to simplify the computation we consider just the first four terms starting with y0 which is set to 1, in this code we identify the Adomian polynomials by A_0, A_1, etc and an image of the code is attached below for easy viewing of Subscripted terms

ClearAll;
Jf[u_, v_, x_] := Gamma[v + 1]/Gamma[u + v + 1] x^(u + v)
Subscript[y, 0] = 1;
Subscript[A, 0] = 1;
Subscript[y, 1][u_, x_] := 2 x + 9/4 Jf[u, 0, x] + Jf[u, 0, x]
Subscript[A, 1][u_, x_] := Subscript[y, 1][u, x]/(
 2 Subscript[y, 0]^0.5)

Subscript[y, 2][u_, x_] := 17/4 Jf[u, 1, x] + (221/32)/Gamma[u + 1] Jf[u, u, x]

Subscript[A, 2][u_, x_] := Subscript[y, 2][u, x]/(2 Subscript[y, 0]^0.5) -  Subscript[y, 1][u, x]^2/(4 Subscript[y, 0]^(3/2))

Subscript[y, 3][u_,x_] := -(9/8) Jf[u, 2, x] + (17^2/(32*Gamma[u + 2]) - 117/(32*Gamma[u + 2])) Jf[u, 1 + u, x] + (3757/(256*Gamma[2 u + ]) - 1521/(512*Gamma[u + 2]^2)) Jf[
    u, 2 u, x]

Subscript[A, 3][u_, x_] := 
 Subscript[y, 3][u, x]/(2 Subscript[y, 0]^0.5) - 
  Subscript[y, 2][u, x]^2/(4 Subscript[y, 0]^(3/2)) - (
  Subscript[y, 1][u, x] Subscript[y, 2][u, x])/(
  2 Subscript[y, 0]^(3/2)) + (3 Subscript[y, 1][u, x]^2)/(
  8 Subscript[y, 0]^(5/2))

Subscript[y, 4][u_, x_] := 
 9/8 Jf[u, 3, x] + 
  1/16 ((-9*34)/(4*Gamma[u + 3]) - (18*17)/(4*Gamma[u + 2]) + (
     27*13)/(4*Gamma[u + 1])) Jf[u, u + 2, x]


Subscript[y, T][u_, x_] := 
 Subscript[y, 0] + Subscript[y, 1][u, x] + Subscript[y, 2][u, x] + 
  Subscript[y, 3][u, x]

Subscript[y, exact][x_] := 
 1/16 E^(-2 x) (E^x + 54 E^(2 x) + 81 E^(3 x) - 
    12 Sqrt[E^(3 x) (1 + 9 E^x)^2])


 Plot[{Subscript[y, T][1.25, x], Subscript[y, T][1.5, x], 
  Subscript[y, T][1.75, x], Subscript[y, T][2, x]}, {x, 0, 1}, 
 PlotTheme -> "Detailed", GridLines -> Automatic, 
 PlotLegends -> {1.25`, 1.5`, 1.75`, 2, Exact}]

Plot[Subscript[y, T][2, x] - Subscript[y, exact][x], {x, 0, 1}, 
 PlotTheme -> "Detailed", GridLines -> Automatic]

Table[{x, Subscript[y, T][2, x] - Subscript[y, exact][x]}, {x, 0, 1, 
   0.1}] // TableForm


 Plot[{Subscript[y, T][1.5, x], poly[x] /. crule // Evaluate}, {x, 0, 
  1}, PlotTheme -> "Detailed", GridLines -> Automatic, 
 PlotLegends -> {Adomian, poly}]