delta[x_, a_] := Block[{}, If[x == a, 1, 0]];
mesh = Table[PDF[NormalDistribution[a, z], x]*
     PDF[NormalDistribution[b, z], y], {a, -1, 1, 0.1},
    {b, -1, 1, 0.1}];
mesh = Flatten[mesh];
cofs = Table[ToExpression[StringJoin[ToString[c],
      ToString[a], ToString[b]]], {a, 0, 20}, {b, 0, 20}];
cofs = Flatten[cofs];
weights = Table[delta[x, a]*delta[y, b], {a, -1, 1, 0.05},
    {b, -1, 1, 0.05}];
weights = Flatten[weights];
expr = cofs . mesh;
residual = D[expr, x] - D[expr, y] - (x*Sin[x*y] -
     y*Sin[x*y]);
integrated = residual*Plus @@ weights;
toMinimize = Plus @@ Flatten[Table[integrated,
      {x, -1, 1, 0.1}, {y, -1, 1, 0.1}]];
toMinimize = toMinimize^2;
cofs = Append[cofs, z];
conds = Abs[(expr /. {x -> 0, y -> 0}) - 1];
opt = NMinimize[{toMinimize + conds}, cofs,
    Method -> {Automatic, PerformanceGoal -> "Speed"}];
u = expr /. Last[opt];
ListPlot3D[Table[u*Plus @@ weights /. {x -> a, y -> b},
   {a, -1, 1, 0.1}, {b, -1, 1, 0.1}]]