How to make a loop with several variables - grounding system?

Hello guys! It's my first time here! I really need your help to solve one doubt! I did write one code to plot some graphs. On first part is OK! It is already done, but, when I try to use one more variable which vary, I have some problems!

This is the code!

close all clear all clc; % Declaração dos parâmetros para simulação % Declaração das constantes elétricas:

    % Declaração da frequência 
    f = logspace(1,7,1000);
    w = 2*pi*f;

    % Permeabilidade Magnética
    mu_0 = 4*pi*1E-7;
    % Permissividade Elétrica do Vácuo
    eps_0 = 8.854*1E-12;
    % Resistividade do Solo (ohm.m)
    rho = 2400;
    % Permissividade Elétrica Relativa
    eps_rel = 15;
    % Permissividade do meio
    eps = eps_rel*eps_0;
    %%%%%%%%%%%%%%%
        % Declaração dos parâmetros do sistema de aterramento:
    % Comprimento do eletrodo em metros
    l = 10:2:100;
    % Raio do eletrodo em metros
    a = (7E-3);
    % Profundidade do eletrodo em metros
    d = 5E-1;
    Fs=20e6; % Frequência de amostragem
    dt=1/Fs; % Intervalo de tempo
    df=500; % Passo de frequências
    ff=0:df:Fs/2; 
    Nf=Fs/df; 
    t=dt*(0:Nf-1); % Força o número de pontos ser o mesmo do número de pontos de frequência

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%DECLARAÇÃO DA FORMA DE ONDA DE CORRENTE %%%        
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        t1=.2e-6;
        t2=40e-6;
        n=2;
        A=t1/t2;
        B=n*t2/t1;
        C=1/n;
        eta=exp(-A*B^C);
        I0=.97e3;

                for aux=1:length(t)
                    A=t(aux)/t1;
                    B=t(aux)/t2;
                    i_t(aux)=(I0/eta)*A^n/(1+A^n)*exp(-B);
                end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

FFT da forma de onda de corrente %%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        
    Iw=fft(i_t);
    Iw=fftshift(Iw);  

    for aux_l = 1:length(l)
    % Cálculo dos parâmetros da linha
    % Resistência do modelo
    resistencia = (rho/pi)*(log((2*l(aux_l))/(sqrt(2*d*a)))-1);
    condutancia = 1/resistencia;
    % Capacitância do modelo
    capacitancia = rho*eps/resistencia; %OK
    % Indutância do modelo
    indutancia = (mu_0/(pi))*(log((2*l(aux_l))/(sqrt(2*d*a)))-1);

    % Impedância Característica
    z_0 = sqrt((j.*w*indutancia)./(condutancia + j.*w*capacitancia));

    % Constante de Propagação
    gama = sqrt(j.*w*indutancia.*(condutancia + j.*w*capacitancia));

    % Impedância    de entrada
    Z = z_0.*coth(gama.*l(aux_l));

    end

The variables is:

'l' is the electrode length... when the variables "resistencia", "condutancia", "capacitancia" and "indutancia" are not variables with the frequency of the signal, that's ok, but, when the variables "rho" and "eps" vary with the frequency, I have tried write a for loop inside a for loop, one to vary frequency and another to vary electrode length.

The frequency variation response follows the relations below and thats ok!

    % Permeabilidade Magnética
    mu_0 = pi*4E-7;
    % Permissividade Elétrica do Vácuo
    eps_0 = 8.854E-12;
    % Resistividade do Solo (ohm.m)
    rho = 300;
    % Condutividade em baixas frequências mS/m  
    sigma_0 = (1/rho)*10^3;
    % Parâmetros constantes (modelo)
    gama_1 = 0.54;
    h = 1.26*sigma_0^(-0.73);

    % Condutividade do meio em função da frequência
    sigma = sigma_0 + sigma_0*h.*(f/1E6).^gama_1;
    % Permissividade Elétrica Relativa
    eps_rel = 12 + (tan(pi*gama_1/2)*(10^-3))/((2*pi*eps_0)*(1E6)^gama_1) .*(sigma_0*h).*f.^(gama_1 - 1);

Thanks for your help!