dt = 0.001;                 % timestep  
tmax = 2;                   % length of stimulus
taumax = 1;                 % maximum timescale of kernel
t=0:dt:tmax;                   % time array as long as stimulus
tau = 0:dt:taumax;               % time-difference array as long as kernel
s=sin(50*t) + 1;            % stimulus can be any function of time

for i = 1:length(tau)       % extend the time vector
    t(end+1) = tmax+dt*i;   % to plot rate as a function of time 
    s(end+1) = 0.0;
end                         % beyond the end of the stimulus

r=zeros(size(t));           % rate: neuron fires as long as stimulus

D=exp(-tau/0.1);            % kernel, D, says how neuron responds to stimulus
imax0 = length(tau);        % number of timepoints in kernel

for i = 1:length(t)         % step through time
    if i < imax0            % number of kernel points to use can not be greater 
        imax = i;           % than number of time points since stimulus started
    else                    % or 
        imax = imax0;       % number of points in the kernel
    end
   
    for j = 1:imax                      % for each point in the kernel
        r(i) = r(i) + D(j)*s(i+1-j);    % sum the kernel multiplied by stimulus that time earlier
    end
    
    r(i) = r(i)*dt;         % Normalization: if dt is small more points are used in the sum, so 
                            % we must scale by this factor so results do
                            % not depend on dt
end

plot(t,r)                   % plot the firing rate.