Unable to plot a solution of NDSolve whilst being able to evaluate it

I want to define a function $$f_c(x,y)$$ as the solution of the PDE $$\Delta f_c(x,y) = g_c(x,y)$$, where $$g_c$$ is a given function depending on a constant $$c$$. I tried the following:

domain = Rectangle[{0, 0}, {1, 1}];
solution[c_] := NDSolve[{Laplacian[h[x, y], {x, y}] == x + c*y, DirichletCondition[h[x, y] == 0, True]}, h, {x, y} \[Element] domain];
mem : f[x_?NumericQ, y_?NumericQ, c_?NumericQ] := mem = Part[Evaluate[h[x, y] /. solution[c]], 1];
f[0.1, 0.2, 3]
Plot3D[f[x, y, 3], {x, y} \[Element] domain]

The fourth line returns the approximate value of $$f_3(0.1,0.2)$$, which to me suggests that everything works. But the plot fails saying NDSolve (...) is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing and a bunch of other similar errors.

I understand the problem probably lies in the way I define $$f$$, but I don't see how to correct it.

This is what ParametricNDSolve is for

domain = Rectangle[{0, 0}, {1, 1}];

Create a parametric function in the parameter c:

pfun = ParametricNDSolveValue[{Laplacian[h[x, y], {x, y}] == x + c*y,
DirichletCondition[h[x, y] == 0, True]},
h, {x, y} \[Element] domain, c];

Use a specific value for c to create the interpolating function:

fun = pfun;
Plot3D[fun[x, y], {x, y} \[Element] domain] • Thanks, but I get problems when I want to take g[x_, y_, c_] := Sum[Part[coefficient, i]*g[x, y, i], {i, 1, c}, with c an integer.. Then I get the error The expression i cannot be used as a part specification when trying to evaluate ParametricNDSolveValue. Nov 18 '21 at 8:26
• How is you comment related to the question? Nov 19 '21 at 6:22
• The function $g_c$ I want to use, uses the parameter $c$ as a summation limit ($x + cy$ was perhaps a wrong MWE). ParametricNDSolveValue does not appear to be able to handle that, hence my question. Is it worthy of a separate post? Nov 19 '21 at 10:51
• @BrazilianCérebro, yes, make new post about that. Nov 19 '21 at 14:30

You are making things way too complicated. I could not understand it. I like to keep things simple myself. Why not just do this?

ClearAll[c, x, y, h, sol];

domain = Rectangle[{0, 0}, {1, 1}];

sol[c_?NumericQ] := NDSolve[{Laplacian[h[x, y], {x, y}] == x + c*y,
DirichletCondition[h[x, y] == 0, True]}, h, Element[{x, y}, domain]];

Plot3D[Evaluate[h[x, y] /. sol], Element[{x, y}, domain]] 