How do I calculate impedance of a (bypass) capacitor (for simulation)?
function [ Z ] = calcCapImp( f,R, L,C )  
Z=sqrt(R^2+(2*pi*f*L-1./(2*pi*C*f)).^2);
end



This is calculating the magnitude of Z, not Z itself.


Z(i)=(Z1*Z2)/(Z1+Z2);%Total impedance is parallel combination



This is combining the two magnitudes as if they are two resistances.

You need to calculate the complex impedance, and combine the complex impedances with the parallel impedances formula

$$Z_{eq} = \frac{Z_1 Z_2}{Z_1 + Z_2}$$

What you calculated was 

$$Z^*_{eq} = \frac{|Z_1| |Z_2|}{|Z_1| + |Z_2|} $$

(* superscript indicates a wrong formula)

which doesn't account for the interaction of phase changes in the two combined impedances, so it doesn't give the correct result.