I need to create a piecewise function that looks like a series of square waves of decreasing width. I will later numerically integrate it, and the code is slow. Is there a better way to do this than with a bunch of heaviside functions?

lam = 1.107e-10; Nzp = 1385; Lf = 28.8346e-3; PHIzp = -(0.1257 - 0.0079*1i)* pi ; Rcs = 30e-6; PHIcs = -(33.5219 - 2.1196*1i) * pi ;

syms rho

Uzp = symfun(0,[rho]);

for m = 0:Nzp

    temp = 2 * (ceil(m/2) - (m/2));
    if temp == 1;
        rmin = ((m - 1/2)*lam*Lf)^(1/2); 
        rmax = ((m + 1/2)*lam*Lf)^(1/2);

        f = symfun(heaviside(rho-rmin),[rho]);
        g = symfun(heaviside(-(rho-rmax)),[rho]);
        h = f * g;

        h2 = h * (exp(1i * PHIzp) - 1);
        Uzp = Uzp + h2;

end

end Uzp = (Uzp + g);

f = symfun(heaviside(-(rho-Rcs)),[rho]); Ucs = f * (exp(1i*PHIcs)-1); Ucs = 1 + Ucs; Uo = Uzp * Ucs;