Why do I get the error message, 'Error using j, Too Many Input/Output Arguments' when I try to execute a function ? I try to combine 21 layer matrice which is 2 by 2. Here is my code. Frequency must change between 300Thz and 500Thz.

function [m11_a,m12_a,m21_a,m22_a] = layer_a(f);

     N_a=2.5;

     delta_za=450e-9;

     k=N_a.*(2*pi.*f/3e8);

     q_a=-k./(2*pi*f*4*pi*10e-7);

        m11_a=cos(k.*delta_za);
        m12_a=((i)./q_a).*sin(k.*delta_za);
        m21_a=(i.*q_a).*sin(k.*delta_za);
        m22_a=cos(k.*delta_za);

     end
----------------------------------------------------------------------------------------------------------------------

     function [m11_b,m12_b,m21_b,m22_b] = layer_b(f);

        N_b=1.5;
        delta_zb=400e-9; 
        k=N_b.*(2*pi.*f/3e8);
        q_b=-k./(2*pi.*f*4*pi*10e-7);

        m11_b=cos(k.*delta_zb);
        m12_b=(i./q_b).*sin(k.*delta_zb);
        m21_b=(i.*q_b).*sin(k.*delta_zb);
        m22_b=cos(k.*delta_zb);

     end
----------------------------------------------------------------------------------------------------------------------
      function [m11,m12,m21,m22]=multiplication(m11_a,m12_a,m21_a,m22_a,m11_b,m12_b,m21_b,m22_b);

m11=m11_a.*m11_b+m12_a.*m21_b;

m12=m11_a.*m12_b+m12_a.*m22_b;

m21=m21_a.*m11_b+m22_a.*m21_b;

m22=m21_a.*m12_b+m22_a.*m22_b;

        end
----------------------------------------------------------------------------------------------------------------------

    function[m11,m12,m21,m22] = ab_multiplication(m11,m12,m21,m22)

[m11,m12,m21,m22]=multiplication(m11,m12,m21,m22,m11,m12,m21,m22);

    end
----------------------------------------------------------------------------------------------------------------------
function [m11_n,m12_n,m21_n,m22_n]= Nemativ_liquid_crystal(f)

        N_n=1.76;
        delta_zn=1e-6; 
        k=N_n.*(2*pi.*f/3e8);
        q_n=-k./(2*pi.*f*4*pi*10e-7);

        m11_n=cos(k.*delta_zn);
        m12_n=(i./q_n).*sin(k.*delta_zn);
        m21_n=(i.*q_n).*sin(k.*delta_zn);
        m22_n=cos(k.*delta_zn);

end ---------------------------------------------------------------------------------------------------------------------- f=3e14:1e12:5e14;

 n=20;

[m11_a,m12_a,m21_a,m22_a] = layer_a(f);

[m11_b,m12_b,m21_b,m22_b] = layer_b(f);

[m11_n,m12_n,m21_n,m22_n]= Nemativ_liquid_crystal(f);

[m11,m12,m21,m22]=multiplication(m11_a,m12_a,m21_a,m22_a,m11_b,m12_b,m21_b,m22_b);

for r=1:5

   [m11,m12,m21,m22] = ab_multiplication(m11,m12,m21,m22);

end

[m11_,m12_,m21_,m22_]=multiplication(m11_b,m12_b,m21_b,m22_b,m11_a,m12_a,m21_a,m22_a);

for r=1:5

   [m11_,m12_,m21_,m22_] = ab_multiplication(m11_,m12_,m21_,m22_);

end

[m11,m12,m21,m22]=multiplication(m11,m12,m21,m22,m11_n,m12_n,m21_n,m22_n);

[m11,m12,m21,m22]=multiplication(m11,m12,m21,m22,m11_,m12_,m21_,m22_);

t(f)=2./((m11+m22)+j(m12-m21));

plot(f,t) ----------------------------------------------------------------------------------------------------------------------

 *My aim is try to combine this matrices which has different refractive indices. And as a result i have to take the graphic which has two parameters transmission coefficients(denoted as t(f) in the code) and frequency. This is too important for my project.
Thank you for your help..

Best regards.. *