K = 1;
Ti = 1;
Td = 1;
Tt = sqrt(Ti * Td);
h = .01;
N = 8;
b = 1;
outputLoLimit = -10.0;
outputHiLimit = 10.0;

bi = K * h / Ti;
ad = ( 2 * Td - N * h) / (2 * Td + N * h);
bd = (2 * K * N * Td) / (2 * Td + N *h);
a0 = h / Tt;

lastDTerm = 0;
ITerm = 0;
ITermOld = 0;

setPoint = 0;

lastY = 0;
Y = [ zeros(1,10) ones(1,90) ]; 
for i = 1:100
    PTerm = K * (b * setPoint - Y(i));
    DTerm = ad *lastDTerm - bd * (Y(i) - lastY);
    
    output1(i) = PTerm * ITerm * DTerm;
    if ( output1(i) > outputHiLimit )
        output2(i) = outputHiLimit;
    elseif (output1(i) < outputLoLimit )
        output2(i) = outputLoLimit;
    else
        output2(i) = output1(i);
    end
    
    ITerm = ITerm + bi*(setPoint - Y(i)) + a0 * (output2(i) - output1(i));
    
    lastDTerm = DTerm;
    lastITerm = ITerm;
end

time = 0:h:99*h;
plot(time, output1, time, output2);