How can I numerically solve for fractional functions and fractional derivatives?

Fractional Iterates

A way to obtain an approximate fractional iterate of a function is to use its Carleman matrix, which is formed from its Taylor coefficients, and then taking the appropriate $p$-th power of the matrix to obtain the series coefficients.

Note that I never said that $p$ had to be an integer; in the example given in the OP, then, we can take the square root of the Carleman matrix to get the coefficients of the semiiterate.

Using the routine CarlemanMatrix[] from this answer, we can do this:

x0 = 0; n = 21; (* expansion point and order *)
sinCM = N[CarlemanMatrix[Sin[x], {x, x0, n}], 30];
shalfCoeffs = MatrixPower[Transpose[sinCM], 1/2, UnitVector[n + 1, 2]];
shalf[x_] = Fold[(#1 x + #2) &, 0, Reverse[shalfCoeffs]];

Notes:

1. I used arbitrary precision to preserve some accuracy, as taking the square root of an unsymmetric matrix can be unstable.
2. I used the "action" form of MatrixPower[], as I only need the coefficients of the semiiterate, which reside in the second row of the matrix square root (hence the Transpose[]).
3. Even with the refinement in point 2, the computation takes a fair bit of time; more so if you're using a larger Carleman matrix.

To check that we now have the Taylor coefficients of the semiiterate, we can use ComposeSeries[]:

ComposeSeries[shalf[x] + O[x]^23, shalf[x] + O[x]^23] - Sin[x] // Chop
   O[x]^23

where we have exploited the fact that shalf[x] consists only of odd powers.

One can then also consider constructing a Padé approximant from the series, like so:

shalfpade[x_] = x PadeApproximant[shalf[x]/x, {x, x0, (n - 1)/2}];

(again exploiting the oddness of the function)

Unfortunately, both the Taylor and Padé approximants have a rather limited range of validity:

Plot[{shalf[x], shalfpade[x]}, {x, -3 π/4, 3 π/4}, 
     PlotStyle -> {RGBColor[7/19, 37/73, 22/31], RGBColor[59/67, 11/18, 1/7]}]

There may be a way to construct better approximants from the Taylor coefficients of the semiiterate, but I haven't investigated them yet.

---

Fractional Derivatives

We can use the general formula of the $\alpha$-th derivative of $x^k$,

$$\frac{\mathrm d^\alpha}{\mathrm dx^\alpha}x^k=\frac{\Gamma(k+1)}{\Gamma(k+1-\alpha)}x^{k-\alpha}$$

and substitute this into the Maclaurin series for the sine. This yields

sinDFrac[x_, a_] = Gamma[a + 1] Sum[(-1)^k Binomial[2 k + 1, a]
                                    x^(2 k + 1 - a)/(2 k + 1)!, {k, 0, ∞}]
   x^(1 - a) Binomial[1, a] Gamma[a + 1]
   HypergeometricPFQ[{1}, {1 - a/2, 3/2 - a/2}, -x^2/4]

that is, a ${}_1 F_2$ hypergeometric function.

In the case $\alpha=\frac12$, we obtain a particularly nice form:

sinDHalf[x_] = FullSimplify[FunctionExpand[sinDFrac[x, 1/2]]]
   Sqrt[2] (Cos[x] FresnelC[Sqrt[2 x/π]] + FresnelS[Sqrt[2 x/π]] Sin[x])

Plot the sine and its semiderivative:

Plot[{Sin[x], sinDHalf[x]}, {x, 0, 2 π}, 
     PlotStyle -> {RGBColor[7/19, 37/73, 22/31], RGBColor[59/67, 11/18, 1/7]}]

Use NonlinearModelFit.

Start with some good guess for what the solution operator might look like, apply that to itself and feed that into NonLinearModelFit:

(* Degree of polynomial approximation *)
k = 5;
(* Points to fit against {{x1, Sin[x1]}, ... } *)
pts = {#, Sin[#]}\[Transpose] &@Range[0, Pi/2, 2/(Pi 2 k)] // N;
(* Our guess for a decent model, c[1]x + c[2]x^2 + ... + c[k] x^k *)
operator[x_] = Plus @@ Array[c@# x^# &, {5}];
(* Apply it to itself *)
doubleOperator = operator[operator[x]]
(* and hopefully find a good fit *)
model = NonlinearModelFit[pts, doubleOperator, Array[c, {5}], x];

GraphicsRow[{
  Plot[
   Abs[Sin[x] - model[x]],
   {x, 0, Pi/2},
   PlotLabel -> "Absolute fit error",
   PlotRange -> All],
  Plot[
   {Sin[x], operator[x] /. model["BestFitParameters"]},
   {x, 0, Pi/2},
   PlotLabel -> Defer[{Sin, Sqrt[Sin]}]]
  }]