l1 norm with fminunc

I am using an optimization function similar to fminunc in matlab. The optimization function takes as input the function and gradient. The gradient is obligatory.

Now given, x = [x1 x2 x3 ...... xn] y = [y1 y2 y3 ...... yn]

I am solving for x and want it to be close to y. The function is something like this:

f(x) = sum{(x1-y1).^2 + (x2-y2).^2 + (x3-y3).^2 + ..... + (xn-yn).^2} + couple of other terms which are also exactly differentiable.

Now I want to add a regularization term to f(x)

i.e. f(x) = f(x) + sum( | x2-x1 | + | x3-x2 | + | x4-x3 | + .... + | x_n-x_n-1 | )

where | | is l1 norm. If I take l2 norm, it is easy to define gradient. I am unable to define the gradient for l1 norm (which I suppose is not differentiable).

I get good results with true point algorithm in fminunc (without defining the gradient), but it is extremely slow.

Can someone suggest how can I define the gradient for the l1 norm?