data = randn(100, 1);            % generate some data
sigma = 0.2;                     % standard deviation of Gaussian kernel
xmin = -4.0;                     % minimum x value for kernel density
xmax = 4.0;                      % maximum x value for kernel density
dx = 0.05*sigma;                 % step size for kernel density
xg = [-4.0*sigma:dx:4.0*sigma];  % x-axis for single Gaussian kernel  
% single Gaussian kernel:
kernel = exp(-0.5*(xg/sigma).^2)/sqrt(2.0*pi)/sigma;
ng = (length(kernel)-1)/2;           % half the length of the Gaussian
x = [xmin:dx:xmax+0.5*dx];       % x-axis for kernel density
kd = zeros(1, length(x));        % vector for kernel density
for i = 1:length(data)           % for every data value ...                
  xd = data(i);
  % index of data value in kernel density vector:
  inx = round((xd-xmin)/dx)+1;
  % start index for Gaussian in kernel density vector:
  k0 = inx-ng;
  % end index for Gaussian in kernel density vector:
  k1 = inx+ng;
  g0 = 1;                        % start index in Gaussian
  g1 = length(kernel);            % end index in Gaussian
  % check whether left side of Gaussian extends below xmin:
  if inx < ng+1
    % adjust start indices accordingly:
    k0 = 1;
    g0 = ng-inx+1;
  end
  % check whether right side of Gaussian extends above xmax:
  if inx > length(kd)-ng
    % adjust end indices accordingly:
    k1 = length(kd);
    g1 = length(kernel)-(inx+ng-length(kd));
  end
  % add Gaussian on kernel density:
  kd(k0:k1) = kd(k0:k1) + kernel(g0:g1);
end
kd /= length(data);               % normalize by number of data points

% plot kernel density:
plot(x, kd)
xlabel('x')
ylabel('Probability density')