Suppose a (noisy) series y, of length m is given. The observations have been sampled at equal distances.
We want to fit a smooth series z to y. We then have to balance two conflicting goals: (1) fidelity to the data and (2) roughness of z. The smoother z is, the more it will deviate from y. Roughness of z can be expressed with differences.

To avoid a lot of tedious algebra, it is advantageous to introduce matrices and vectors
$$ Q=|\mathbf{y}-\mathbf{z}|^{2}+\lambda|\mathbf{D z}|^{2}$$ (1)
where the first term $$|\mathbf{y}-\mathbf{z}|^{2}$$ indicates the quadratic norm of error, and $\mathbf{D} $ is a matrix such that $\mathbf{D} \mathbf{z}=\Delta \mathbf{z}$.

To solve the z, such that the $$ Q=|\mathbf{y}-\mathbf{z}|^{2}+\lambda|\mathbf{D z}|^{2}$$ can be minimized, we should introduce the partial derivatives( the minimal or maximal is always at the point when the partial derivatives be zero).
$$ \partial Q / \partial \mathbf{z}^{\prime}=-2(\mathbf{y}-\mathbf{z})+2 \lambda \mathbf{D}^{\prime} \mathbf{D} \mathbf{z} $$
and equating this to 0 we get the linear system of equations
$$ \left(\mathbf{I}+\lambda \mathbf{D}^{\prime} \mathbf{D}\right) \mathbf{z}=\mathbf{y} $$ where I is the identity matrix.

The matrix $\mathbf{D} $ is actually an differential matrix. To make more sense, In the avbove, For example,when m=5, $\mathbf{D} $ would be
$$ \mathbf{D}=\left[ \begin{array}{rrrrr}{-1} & {1} & {0} & {0} & {0} ;\ {0} & {-1} & {1} & {0} & {0} ;\ {0} & {0} & {-1} & {1} & {0}; \ {0} & {0} & {0} & {-1} & {1}\end{array}\right]$$

Implementation in Matlab is very straightforward, because we can use the built-in function diff() to compute $D$:

m = length(y);
E = eye(m);
D = diff(E);
z =  (E + lambda * D' * D)\y;

A better implementation can be download at https://www.mathworks.com/matlabcentral/fileexchange/25634-smoothn

 %--- Example #1: smooth a curve ---
  x = linspace(0,100,2^8);
  y = cos(x/10)+(x/50).^2 + randn(size(x))/10;
  y([70 75 80]) = [5.5 5 6];
  z = smoothn(y); % Regular smoothing
  zr = smoothn(y,'robust'); % Robust smoothing
  subplot(121), plot(x,y,'r.',x,z,'k','LineWidth',2)
  axis square, title('Regular smoothing')
  subplot(122), plot(x,y,'r.',x,zr,'k','LineWidth',2)
  axis square, title('Robust smoothing')

Result of Example1

Reach matcodelab@gmail.com