Specify non-homogeneous initial conditions - Partial Differential Equation Toolbox

Hi,

I am using the Partial Differential Equation Toolbox to solve one PDE over a two-dimensional quadratic domain with a mesh composed of 3-by-3 nodes (for the sake of simplicity).

I would like to assign the following initial conditions: init = [2 3 4; 5 6 7; 8 9 10] - meaning that the value of 2 will be assigned to the node located in the upper left corner of my geometry, and so on.

I am using the command setInitialConditions(model,init) to set up the initial conditions but I do not know how to pass the init matrix correctly. If I simply pass it as reported in the command above Matlab rejectes it. I think I would need to create a local function to pass the matrix init, but there is no clear example in the documentation page for setInitialConditions https://se.mathworks.com/help/pde/ug/pde.pdemodel.setinitialconditions.html

I am trying to pass the initial conditions as follow

% Initial conditions
setInitialConditions(model,@initfun);

function g = initfun(location)
   nr = length(location.x); % Number of nodes
   g = zeros(1,nr); % Allocate f
   init = [2 3 4; 5 6 7; 8 9 10];
   g(1,:) = init(:);
end

But I get the following error

Error using
pde.GeometricInitialConditions/checkIcFcnHdlArgCounts (line
230)
Function handle specifying an initial condition must accept
one input argument and return one output argument.

Error in pde.GeometricInitialConditions/precheckIC (line 207)
         pde.GeometricInitialConditions.checkIcFcnHdlArgCounts(icval);

Error in pde.GeometricInitialConditions/set.InitialValue
(line 118)
      self.precheckIC(val);

Error in pde.GeometricInitialConditions (line 75)
          obj.InitialValue = varargin{2};

Error in pde.PDEModel/setInitialConditions (line 136)
            ic =
            pde.GeometricInitialConditions(argsToPass{:});

Error in feature_testing (line 45)
setInitialConditions(model,@initfun);

Additionally, if I modify the function in the following form

function g = initfun(location)
   nr = length(location.x); % Number of nodes
   g = zeros(1,nr); % Allocate f
   location.x
   g(1,:) = ones(1,length(location.x));
end

so to have the output of location.x I get the following output

ans =

     0


ans =

     0


ans =

     0


ans =

     0


ans =

     0


ans =

     0


ans =

  Columns 1 through 10

     1     1     1     1     2     2     2     2     3     3

  Columns 11 through 16

     3     3     4     4     4     4

>> 

I do not understand where the six zeros are coming from since the x location of my nodes is represented by the 1 1 1... 4 4 4 part.

Ideas?

Many thanks,

Jacopo