## backward finite difference with fractional differential equat. (G-L definition)

Hi everyone!
I'm having a problem with backward finite difference method. I'm dealing with this problem using Mathematica software (and I'm going to paste some code below, maybe you don't exactly know the syntax, but it's pretty easy I think)..

I need to use the backward finite difference method with the Grunwald-Letnikov fractional derivative definition (probably it's an uncommon topic..).
Long story short, G-L definition for the fractional derivative is this(two links here):

http://en.wikipedia.org/wiki/Gr%C3%B...kov_derivative
and
http://www.xuru.org/fc/GrunwaldLetnikov.asp

and can be usually used for backw.FD in order to approximate some fractional differential equation.

This is my problem.
I have a first order mass-spring differetial equation and and I have to write a code in Mathematica to approximate it (with backwFD method) and then compare it with the fractional differential equation case. (ps: the same problem appears in the second order equation (mass-spring-damper) approximation, but let's just limit this on the first order).

The problem is that, for the ordinary case I have no problem in using FD method, instead I have a big problem in using FD method for the fractional case: I'll write some code to be more precise.

---
first order - ordinary
---
Code:
massa = 0;cvisc = 5000;kappa = 1000;forzante[t_] := 2000*UnitStep[t];
LL = 30;
nn = 31;
deltah = LL/(nn - 1)
ii = Range[0, nn];
tt = deltah*(ii - 1);
Clear[u]; Clear[eq]; Clear[eqi]; Clear[ui];
ui = Array[u, nn + 1, 0];
eqi = Array[eq, nn + 1, 0];

For[i = 1, i <= nn, i++,
eq[i] = cvisc*((u[i] - u[i - 1])/(deltah)) + kappa*(u[i]) ==
forzante[tt[[i + 1]]]]
eq[0] = u[1] == 0;

solApprox := NSolve[eqi, ui]; uisol =
Table[ui[[i]] /. solApprox[[1]], {i, nn + 1}];
dataApprox = Table[{tt[[i]], uisol[[i]]}, {i, 2, nn + 1}];
plot1 = ListPlot[dataApprox, PlotStyle -> {Purple},
PlotRange -> Automatic]

---
first order - fractional
---
everything as before, but with this FD cicle

Code:
For[i = 1, i <= nn, i++,
eq[i] = (1/deltah)*cvisc*
Sum[(-1)^r*Binomial[.5, r]*u[i - r], {r, 0, i}] + kappa*(u[i]) ==
forzante[tt[[i + 1]]]]
which is made by substituting the first order derivative with a fractional derivative of order 1/2 using GL definition.

----

The first code (ordinary case) works flawlessly, but the second doesn't: the curve is completely different depending on the "deltah" imposed (it gives no problems when deltah=1 , but smaller values are generally used in approximations...and that's where this code gives problems).

Now, i've tried quite everything with this code but to no avail...the problem is for sure given by deltah but i cannot find a way.

It would be awesome if you could give me some hints on this. I think you could give me some hints also if you don't know Mathematica.

PS: Since it's my first post here (by the way...this forum is awesome and I say Hi to everybody :) ), please excuse me if I'm posting it in the wrong section or if I'm going against some netiquette.