TLM code not evolving as needed

Hi All,

I am using MATLAB to code a 2D TLM problem. In essence, it is a mesh grid which I excite at one point and measure at some other point the response. There is also boundary conditions, but most importantly are the scattering and connection steps. My code is not performing as it should. Here is a picture to illustrate the concept:

http://upload.wikimedia.org/wikipedia/commons/thumb/f/f8/SingleNode2DTLM.png/500px-SingleNode2DTLM.png

Finally, I do a discrete Fourier transform and plot versus frequency. The code is not evolving as it should, since V1r or V1i at [let's say] (2,27,27) is 0. Does anyone see what is wrong with the below code? The first for loop is for time, the second and third is for scattering and fourth and fifth for connection. Can I group the 2nd and 3rd with 4th and 5th or do the cells need to populated first? How come the impulses are not travelling to other cells? Anyone have any tips???

clc clear all

a = 0.2;

d = 0.2;

source_j = 26; source_k = 26;

probe_j = 29; probe_k = 45;

L = (4*pi)*10^(-7); C = 0.5*(8.854*10^(-12)); del_L = 0.0002; vL = 1/sqrt(L*C); ZL = sqrt(L/C); del_t = del_L/(2*vL); na = (a/del_L)+1; nd = (d/del_L)+1; nt = 1013;

Vy = zeros (nt, nd, na);

V1i = zeros(nt,nd,na); V2i = zeros(nt,nd,na); V3i = zeros(nt,nd,na); V4i = zeros(nt,nd,na);

V1r = zeros(nt,nd,na); V2r = zeros(nt,nd,na); V3r = zeros(nt,nd,na); V4r = zeros(nt,nd,na); f = linspace(1*(10^9),100*(10^9),100);

F=zeros(1,100);

%excite just one cell
V1i(1,26, 26) = 1;
V2i(1,26, 26) = 1;
V3i(1,26, 26) = 1;
V4i(1,26, 26) = 1;

for i=1:nt %nt-1

    for j = 1:na-1

        for k = 1:nd-1

            V1r(i,j,k) =0.5*( -V1i(i,j,k) + V2i(i,j,k) + V3i(i,j,k) + V4i(i,j,k));
            V2r(i,j,k) = 0.5*(V1i(i,j,k) - V2i(i,j,k) + V3i(i,j,k) + V4i(i,j,k));
            V3r(i,j,k) = 0.5*(V1i(i,j,k) + V2i(i,j,k) - V3i(i,j,k) + V4i(i,j,k));
            V4r(i,j,k) = 0.5*(V1i(i,j,k) + V2i(i,j,k) + V3i(i,j,k) - V4i(i,j,k));

        end
    end

        for j=2:na-1
            for k = 2:nd-1

            V1i(i+1,j,k) = V3r(i,j+1,k);
            V2i(i+1,j,k) = V4r(i,j,k-1);
            V3i(i+1,j,k) = V1r(i,j-1,k);
            V4i(i+1,j,k) = V2r(i,j,k+1);
        end
    end

    for k =1 % zmin
        V2i(i+1,j,k) = -V2r(i,j,k);
    end

    for j = 1
        V3i(i+1,j,k) = -V3r(i,j,k);
    end

            for k = nd 
                V4i(i+1,j,k) = -V4r(i,j,k);
            end

            for j = na
                V1i(i+1,j,k) = -V1r(i,j,k);
            end

    Vy(i,probe_j,probe_k) = 0.5*(V1i(i,probe_j,probe_k) + V2i(i,probe_j,probe_k) + V3i(i,probe_j,probe_k) + V4i(i,probe_j,probe_k));

    F=F+Vy(i,probe_j,probe_k)*exp(-2*pi*f);

    V1i(i+1,:,:)=V1i(i,:,:);
    V2i(i+1,:,:)=V2i(i,:,:);       %Denys
    V3i(i+1,:,:)=V3i(i,:,:);
    V4i(i+1,:,:)=V4i(i,:,:);
end            

    plot(f, abs(F));

[Theophanes Raptis]

The posted code does not even work properly.Last 4 loops do nat have correct limits. All internal loops could also be vectorized.