modified
fxy = diffmat[[All, 4, All]]; (* to be multiplied by gridpoint values*) derivxy=MapThread[Join[#1, {#2}] &, {pts, fxy . ijy[[All, -1]]}] (* {{i,j,fxy},...}*) ip = Interpolation[derivxy, InterpolationOrder -> 1] Show[{Plot3D[ip[i, j], Element[{i, j}, mesh],AxesLabel -> {i, j, "fxy"}], Mesh -> {All, Full}[[1]]], Graphics3D[Point[derivxy]]}] 
Hope
Hope it helps!
fxy = diffmat[[All, 4, All]]; (* to be multiplied by gridpoint values*) derivxy=MapThread[Join[#1, {#2}] &, {pts, fxy . ijy[[All, -1]]}] (* {{i,j,fxy},...}*) ip = Interpolation[derivxy, InterpolationOrder -> 1] Plot3D[ip[i, j], Element[{i, j}, mesh],AxesLabel -> {i,j, "fxy"}]
Hope it helps!
modified
fxy = diffmat[[All, 4, All]]; (* to be multiplied by gridpoint values*) derivxy=MapThread[Join[#1, {#2}] &, {pts, fxy . ijy[[All, -1]]}] (* {{i,j,fxy},...}*) ip = Interpolation[derivxy, InterpolationOrder -> 1] Show[{Plot3D[ip[i, j], Element[{i, j}, mesh],AxesLabel -> {i, j, "fxy"}, Mesh -> {All, Full}[[1]]], Graphics3D[Point[derivxy]]}]
Hope it helps!
pts = ijy[[All, 1 ;; 2]]; mesh = DiscretizeRegion[ResourceFunction["NonConvexHullMesh"][pts,2]]2]]; Show[mesh, Graphics[Point[pts]]] 
fxy = diffmat[[All, 4, All]]; (* to be multiplied by gridpoint values*) derivxy=MapThread[Join[#1, {#2}] &, {pts, fxy . ijy[[All, -1]]}] (* {{i,j,fxy},...}*) ListPlot3D[derivxyip = Interpolation[derivxy, MeshInterpolationOrder -> All1] Plot3D[ip[i, PlotRangej], ->Element[{i, Allj}, mesh],AxesLabel -> {i, j, "fxy"}] & 
attention: ListPlote3D incorrectly shows interpolation using ConvexHullMesh (not the mesh mesh we defined)
pts = ijy[[All, 1 ;; 2]]; mesh = DiscretizeRegion[ResourceFunction["NonConvexHullMesh"][pts,2]]
fxy = diffmat[[All, 4, All]]; (* to be multiplied by gridpoint values*) derivxy=MapThread[Join[#1, {#2}] &, {pts, fxy . ijy[[All, -1]]}] (* {{i,j,fxy},...}*) ListPlot3D[derivxy, Mesh -> All, PlotRange -> All,AxesLabel -> {i, j, "fxy"}] &
attention: ListPlote3D incorrectly shows interpolation using ConvexHullMesh (not the mesh mesh we defined)
pts = ijy[[All, 1 ;; 2]]; mesh = DiscretizeRegion[ResourceFunction["NonConvexHullMesh"][pts,2]]; Show[mesh, Graphics[Point[pts]]]
fxy = diffmat[[All, 4, All]]; (* to be multiplied by gridpoint values*) derivxy=MapThread[Join[#1, {#2}] &, {pts, fxy . ijy[[All, -1]]}] (* {{i,j,fxy},...}*) ip = Interpolation[derivxy, InterpolationOrder -> 1] Plot3D[ip[i, j], Element[{i, j}, mesh],AxesLabel -> {i,j, "fxy"}]
code (we need a nonconvex mesh for your data)
pts = ijy[[All, 1 ;; 2]]; mesh = DelaunayMesh[ pts];DiscretizeRegion[ResourceFunction["NonConvexHullMesh"][pts,2]]
nf = Nearest[pts -> "Index"]; neighbors = Table[nf[pts[[ii]], 8 + 1], {ii, 1, Length[pts]}]; diffmat = Table[index = neighbors[[i]]; {xi, yi} = pts[[index[[1]]]];(*aktueller Punkt*) nj = Length[index] - 1;(*aktuelle Anzahl Nachbarn*) sparse = SparseArray[ Flatten[Table[{{j, index[[j + 1]]} -> 1, {j, index[[1]]} -> -1}, {j, 1, nj}], 1], {nj, Length[pts]}]; pinv = PseudoInverse[ Table[{xj, yj} = pts[[index[[j + 1]]]]; {xj - xi, yj - yi, 1/2 (xj - xi)^2, (xj - xi) (yj - yi), 1/2 (yj - yi)^2, 1/6 (xj - xi)^3, 1/2 (xj - xi)^2 (yj - yi), 1/2 (xj - xi) (yj - yi)^2, 1/6 (yj - yi)^3}, {j, 1, nj}]]; pinv . sparse, {i, 1, Length[pts]}]; fxy = diffmat[[All, 4, All]]; (* to be multiplied by gridpoint values*) MapThread[Join[#1derivxy=MapThread[Join[#1, {#2}] &, {pts, fxy . ijy[[All, -1]]}] // (* {{i,j,fxy},...}*) ListPlot3D[# ListPlot3D[derivxy, Mesh -> All, PlotRange -> All,AxesLabel -> {i, j, "fxy"}] & 
attention: ListPlote3D incorrectly shows interpolation using ConvexHullMesh (not the mesh mesh we defined)
code
pts = ijy[[All, 1 ;; 2]]; mesh = DelaunayMesh[ pts]; nf = Nearest[pts -> "Index"]; neighbors = Table[nf[pts[[ii]], 8 + 1], {ii, 1, Length[pts]}]; diffmat = Table[index = neighbors[[i]]; {xi, yi} = pts[[index[[1]]]];(*aktueller Punkt*) nj = Length[index] - 1;(*aktuelle Anzahl Nachbarn*) sparse = SparseArray[ Flatten[Table[{{j, index[[j + 1]]} -> 1, {j, index[[1]]} -> -1}, {j, 1, nj}], 1], {nj, Length[pts]}]; pinv = PseudoInverse[ Table[{xj, yj} = pts[[index[[j + 1]]]]; {xj - xi, yj - yi, 1/2 (xj - xi)^2, (xj - xi) (yj - yi), 1/2 (yj - yi)^2, 1/6 (xj - xi)^3, 1/2 (xj - xi)^2 (yj - yi), 1/2 (xj - xi) (yj - yi)^2, 1/6 (yj - yi)^3}, {j, 1, nj}]]; pinv . sparse, {i, 1, Length[pts]}]; fxy = diffmat[[All, 4, All]]; (* to be multiplied by gridpoint values*) MapThread[Join[#1, {#2}] &, {pts, fxy . ijy[[All, -1]]}] // ListPlot3D[#, Mesh -> All, PlotRange -> All,AxesLabel -> {i, j, "fxy"}] &
code (we need a nonconvex mesh for your data)
pts = ijy[[All, 1 ;; 2]]; mesh = DiscretizeRegion[ResourceFunction["NonConvexHullMesh"][pts,2]]
nf = Nearest[pts -> "Index"]; neighbors = Table[nf[pts[[ii]], 8 + 1], {ii, 1, Length[pts]}]; diffmat = Table[index = neighbors[[i]]; {xi, yi} = pts[[index[[1]]]];(*aktueller Punkt*) nj = Length[index] - 1;(*aktuelle Anzahl Nachbarn*) sparse = SparseArray[ Flatten[Table[{{j, index[[j + 1]]} -> 1, {j, index[[1]]} -> -1}, {j, 1, nj}], 1], {nj, Length[pts]}]; pinv = PseudoInverse[ Table[{xj, yj} = pts[[index[[j + 1]]]]; {xj - xi, yj - yi, 1/2 (xj - xi)^2, (xj - xi) (yj - yi), 1/2 (yj - yi)^2, 1/6 (xj - xi)^3, 1/2 (xj - xi)^2 (yj - yi), 1/2 (xj - xi) (yj - yi)^2, 1/6 (yj - yi)^3}, {j, 1, nj}]]; pinv . sparse, {i, 1, Length[pts]}]; fxy = diffmat[[All, 4, All]]; (* to be multiplied by gridpoint values*) derivxy=MapThread[Join[#1, {#2}] &, {pts, fxy . ijy[[All, -1]]}] (* {{i,j,fxy},...}*) ListPlot3D[derivxy, Mesh -> All, PlotRange -> All,AxesLabel -> {i, j, "fxy"}] &
attention: ListPlote3D incorrectly shows interpolation using ConvexHullMesh (not the mesh mesh we defined)