Simulating PID controller in Matlab

I am trying to explore PID control with a simulation. My problem is that the derivative term does not behave as expected in my simulation: No matter how large its coefficient it seems unable to prevent overshoots.

I tried to tune the system with the Ziegler-Nichols method and using the parameters the should produce "no overshoot" according to the Wikipedia entry of the Ziegler-Nichols method.

My "process" is simply the 5-sample exponential moving average of the PID controller output. Controller output is limited to [-1,1]

Can anybody see what the problem with this code is:

PVset = 0.5;            % set point for process variable

steps = 100;          % time periods/steps to simulate

PV = 0;               % process variable (time series)
error = 0;            % error in PV (time series)
der = 0;              % derivative of error (time series)

Ku = 1;
Tu = 2;

Kp = Ku*.2;
Ki = 2*Kp/Tu;
Kd = Kp*Tu/3;

for k=2:steps
    % calculate error
    error(end+1) = PVset-PV(end);
    % P term
    delta = Kp*error(end);
    % I term
    delta = Ki*sum(error) + delta;
    % D term
        % instantaneous derivative time series
    der = error(2:end)-error(1:end-1);
        % average derivative over x periods
    per = 4;
    if k >= per+1
        avgDer = mean(der(end-(per-1):end));
        delta = Kd*avgDer+delta;
    end
    % limit PID controller output to -1 -> 1
    if abs(delta) > 1
        delta = sign(delta);
    end
    % simulate process behaviour
    p = 5;
    PV(k+1) = (p-1)/p*PV(k-1) + 1/p*delta;
end

figure
hold on
plot((0:steps), PV)
plot((1:steps), error, 'r')