How to transfer a solution of a PDE as the coeeficient for the next PDE?

Hi everyone,

I need to solve a system of PDEs for an optimization problem. The system I would need to solve is the original PDE and the related Adjoint equation in these specific forms over the unit square in 2D.

PDE: \nabla * ( c \nabla u) + a*u^3= f (Solve for u)

ADJOINT \nabla * ( c \nabla p) + a* u^2* p = u (Solve for p)

c,a and f are functions.

My Code (so far) is

%Generate PDE
model=createpde();
g=geometrysquare();
geometryFromEdges(model,g);
generateMesh(model,'Hmax',0.01);

%Generate Conditions and Coefficients
c=@(location,state)noise(0.5,20,location.x,location.y);
a=@(location,state)noise(0.5,20,location.x,location.y).*state.u.^2;
f1=@(location,state) sin(pi.*location.x).*sin(pi.*location.y).^2;
applyBoundaryCondition(model,'dirichlet','Edge',1:model.Geometry.NumEdges,'u',0);
specifyCoefficients(model,'m',0,'d',0,'c',c,'a',a,'f',f1);

%Solve first PDE
results=solvepde(model);
pdeplot(model,'XYData',results.NodalSolution) %just to see how it looks
axis equal

%Generate a function which is evaluated in mesh point 
v=linspace(0,1,100); %Create Linspace with mesh size related to generateMesh
[X,Y] = meshgrid(v);
querypoints = [X(:),Y(:)]';
A = interpolateSolution(results,querypoints);
b=length(X);
A = reshape(A,size(X));
A(1,:)=zeros(b,1); %Since MATLAB is not perfectly exact, we will make the edges 0 by default
A(b,:)=zeros(b,1);
A(:,1)=zeros(1,b);
A(:,b)=zeros(1,b);

%create adjoint Model
model2=createpde();
geometryFromEdges(model2,g);
generateMesh(model2,'Hmax',0.1);

%Create adjoint coefficient functions
f2=@(location,state) A(sign(location.x)*ceil(location.x.*b),sign(location.y)*ceil(location.y.*b)); 
a2=@(location,state) noise(0.5,20,location.x,location.y).*A(sign(location.x).*ceil(location.x.*b),sign(location.y).*ceil(location.y.*b)).^2;


% Put this function as coefficient f for the adjoint equation
applyBoundaryCondition(model2,'dirichlet','Edge',1:model2.Geometry.NumEdges,'u',0);
CB=specifyCoefficients(model2,'m',0,'d',0,'c',1,'a',a,'f',f2);

% Solve the adjoint equation
resultsadjoint=solvepde(model2);

% Save the adjoint solution similar to the original solution at the same points
B=interpolateSolution(resultsadjoint,querypoints);
B=reshape(B,size(X));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%The function creates geometry for the square (0,1)^2 to be inserted for the PDE later
function g=geometrysquare()

% Coordinates
lowerLeft  = [0  ,0  ];
lowerRight = [ 1 , 0  ];
upperRight = [1 , 1];
upperLeft =  [0 , 1];
% Geometry matrix
S = [3,4 lowerLeft(1), lowerRight(1), upperRight(1), upperLeft(1), ...
         lowerLeft(2), lowerRight(2), upperRight(2), upperLeft(2)];                     
gdm = S';
% Names
ns = 'S';
% Set formula 
sf = 'S';
% Invoke decsg
g = decsg(gdm,ns,sf');
end

Then I get an error message "Index in position 1 is invalid. Array indices must be positive integers or logical values." while evaulating my matrix A in location.x . But why does this error occur?

Normally ceil(x) rounds UP to the next interger. Therefore, my smallest index should be 1 and not 0.

Hope that someone can help me here.

All the Best!