freq : 0.05;
stepSize : 2e-2;
stepLim : ceiling(1/(freq*stepSize));
hyst : 0.15;
Xgrid : [0];
Xaxis : makelist(k*stepSize,k,0,stepLim)$
refSignal(t) := sin(2*%pi*freq*t);
intSignal : [0];
errVal : [0];
  /* this loop takes quite some time to execute */
pwmSignal : block(x:[1],for i:1 thru stepLim step 1 do block(
    errVal : append(errVal,[refSignal(i*stepSize) - x[i]]),
    intSignal : append(intSignal,[intSignal[i]+(errVal[i]*stepSize)]),
    if (intSignal[i+1] > hyst) then 
        x : append(x, [1])
    elseif(intSignal[i+1] < (-hyst)) then 
        x : append(x, [-1])
    else
        x : append(x, [x[i]])),
    x
)$
Xgrid : block(x:[],for k:2 thru stepLim step 1 do 
    if(pwmSignal[k-1]#pwmSignal[k]) then 
        x : append(x,[Xaxis[k]]),
    setify(x));
plot1 : gr2d(
    user_preamble="set format x \"\";set bmargin 0;
        set format y \"\"",
    font="Times",grid=true,xtics=Xgrid,
    points_joined=true,point_type=dot,
    points(Xaxis,errVal),yrange=[-2,2],
    color=green,explicit(refSignal(t),t,0,20),
    ytics={-2,-1,0,1,2},ylabel="Ref. and Error"
)$
plot2 : gr2d(
    user_preamble="set format x \"\"; set tmargin 0;
        set bmargin 0",
    font="Times",grid=true,xtics=Xgrid,
    points_joined=true,point_type=dot,
    color=magenta,points(Xaxis,intSignal),
    yrange=[-.25,.25],ytics={-.2,0,.2},
    ylabel="Integration"
)$
plot3 : gr2d(
    user_preamble="set format x \"\";set bmargin 1.75;
        set xlabel \"Time\" 0,1.5",
    font="Times",grid=true,xtics=Xgrid,
    points_joined=true,point_type=dot,
    points(Xaxis,pwmSignal),yrange=[-1.25,1.25],
    ytics={-1,0,1},ylabel="Sigma-Delta PWM"
)$
draw(terminal=svg,file_name="Sigma-delta PWM",
    plot1,plot2,plot3);