Arrayfun question - getting error

Hello, I have a question with regard to using the arrayfun function. I am trying to perform a simple stencil computation using a nested function as shown below.

However, I am getting an error:

Attempted to access u(1,2,2); index out of bounds because size(u)=[1,10,1].

Error in lappy/loop (line 23)
            n = (u(i-1,j,k) + u(i+1,j,k) + u(i,j-1,k) + u(i,j+1,k) + u(i,j,k-1) +
            u(i,j,k+1))*(1/6);

Error in lappy (line 29)
    un = arrayfun(@loop,i,j,k);

I have attached the code below.

function [u,u1] = lappy()
  nx = 10;
  ny = 10;
  nz = 10;
  iters = 10;

    u = zeros(nx,ny,nz);
    u(1,:,:) = 1.0;
    u(nx,:,:) = 1.0;
    u(:,1,:) = 1.0;
    u(:,ny,:) = 1.0;
    u(:,:,1) = 1.0;
    u(:,:,nz) = 1.0;

    i = 1:nx;
    j = 1:ny;
    k = 1:nz;

      function [n] = loop(i,j,k)
          if (i == 1 || j == 1 || k == 1 || i == nx || j == ny || k == nz)
              n = u(i,j,k);
          else
              n = (u(i-1,j,k) + u(i+1,j,k) + u(i,j-1,k) + u(i,j+1,k) + u(i,j,k-1) + u(i,j,k+1))*(1/6);
          end
      end

    tic;
    for t = 1:iters
        un = arrayfun(@loop,i,j,k);
        u = un;
    end
    toc;

end

Answer 1

[Guillaume]

The error you're getting has nothing to do with arrayfun. As it tells you the size of u is [1,10,1], so it looks like the code you've posted is not exactly the code you're running (since it looks like nx and nz are both 1).

In any case, your arrayfun is only going to work if i, j and k are all the same size (otherwise it will throw an error) and if it is, is not going to work the way you want. On the first pass, arrayfun will pass the first elements of i, j and k to your loop function. On the second pass, it will pass the second elements of i, j and k, etc., on the nth pass it passes the nth elements of all the arrays. That is arrayfun moves through all the arrays at the same time. It never passes the 1st element of the 1st array with the 2nd element of the second.

So if you want to test all the combination of i, j, k, you either have to write a triple for loop, or use ndgrid to generate these combinations:

[i, j, k] = ndgrid(1:nx, 1:ny, 1:nz); %and i,j,k are poor variable names
un = arrayfun(@loop, i, j, k);