Summary:
There are infinitely many functions which satisfy $\left(1+\partial_x\partial_y\right)f(x,y)=0$, that is why Mathematica does not immediately give you the specific answer $sin(x+y)$. For example both $$-J_0\left(2 \sqrt{x} \sqrt{y}\right)\quad,\quad-\frac{2 \sqrt{y} J_1\left(2 \sqrt{x} \sqrt{y}\right)}{\sqrt{x}}$$ are functions that satisfy given differential equation.
Code to solve homogeneous differential equations symbolically
Following code finds a series expansion for any homogeneous differential equation:
ClearAll[symbolicDSolve];
symbolicDSolve = Function[{equation, function},
Module[{variables, exponent, exponents, sum, sum2, c, n, a,
equation2, recursion, result, boundary, finalequations,
finalresult},
variables = List @@ function;
exponents = Table[Unique["n"], Length[variables]];
exponent =
If[Length[exponents] == 1, Last[exponents], exponents];
sum[a_ + b_] := sum[a] + sum[b];
sum[x_^(n_ + a_)] :=
sum[x^n] /; And[MemberQ[variables, x], MemberQ[exponents, n]];
sum[b_ x_^(n_ + a_)] :=
sum[(b /. n -> n - a) x^n] /;
And[MemberQ[variables, x], MemberQ[exponents, n]];
simplify = {sum2[a_ ] + sum2[b_ ] :> sum2[a + b],
a_ sum2[b_] :> sum2[a b]};
equation2 = (equation //.
function :>
sum[(c @@ exponents) (Times @@ (variables^exponents))]) //.
Evaluate[Derivative[a__][Head[function]][##] & @@ variables] :>
sum[(D[(c @@
exponents) (Times @@ (variables^exponents)), ##] & @@
Transpose[{variables, {a}}])];
recursion =
Simplify[
Coefficient[
Last[equation2 /. sum -> sum2 //.
simplify], (Times @@ (variables^exponents))]];
result =
Evaluate[(Inactive[
Sum][(c @@
exponents) (Times @@ (variables^exponents)), ##] & @@
Table[{exponents[[i]], 0, Infinity}, {i,
Length[variables]}])] /.
Flatten[{RSolve[recursion == 0, (c @@ exponents), exponent]}]
]
];
Analysis of the original question:
Your question:
In[5]:= symbolicDSolve[
D[f[x, y], x, y] + f[x, y], f[x, y]]
Out[5]= Inactive[Sum][((-1)^(-1 + n13) x^n13 y^n14 C[1][n13 - n14])/(
Pochhammer[2, -1 + n13] Pochhammer[2 - n13 + n14, -1 + n13]), {n13,
0, \[Infinity]}, {n14, 0, \[Infinity]}]
which is simply
$$\underset{\text{n13}=0}{\overset{\inf }{\sum }}\underset{\text{n14}=0}{\overset{\inf}{\sum }}\frac{(-1)^{\text{n13}-1} x^{\text{n13}}y^{\text{n14}}c_1(\text{n13}-\text{n14})}{(2)_{\text{n13}-1}(-\text{n13}+\text{n14}+2)_{\text{n13}-1}}$$
for undetermined coefficients $c(n13-n14)$. We see that differential equation as a constraint reduced the two dimensional unknown $c(n13,n14)$ to one dimensional unknown $c(n13-14)$ ; still, there are infinitely many solutions!
Let us assume that $c(0)=1$ and all other $c(i)=0$. That means, out solution is $$f(x,y)=\underset{\text{a}=0}{\overset{\infty }{\sum }}\frac{(-1)^{a+1} c_1(0) x^a y^a}{\Gamma (a+1)^2}$$
We can let mathematica calculate this for us, we immediately get the Bessel function:
In[8]:=Sum[((-1)^(1 + a)x^ay^aGamma[2]C[1][0])/(Gamma[1+a]Gamma[1+a]),
{a,0,Infinity}]Out[8]= -BesselJ[0, 2 Sqrt[x] Sqrt[y]] C[1][0]
Of course we can directly check that the result satisfies original differential equation:
In[709]:= (Derivative[1, 1][f][x, y] + f[x, y]) /.
f -> Function[{x, y}, -((2 Sqrt[y] BesselJ[1, 2 Sqrt[x] Sqrt[y]])/
Sqrt[x])] // FullSimplify
Out[709]= 0
Other examples of the code:
First order, one variable: $a f'(x)+b f(x)=0$:
In[16]:= symbolicDSolve[a D[f[x], x] + b f[x], f[x]] // Activate
Out[16]= -((a E^(-((b x)/a)) C[1])/b)
$$f(x)=-\frac{a c_1 e^{-\frac{b x}{a}}}{b}$$
Second order, one variable: $f''(x)+a x f'(x)+ b f(x)=0$:
In[26]:=symbolicDSolve[f''[x]+axf'[x]+bf[x],f[x]]
//Activate//FullSimplify
Out[26]= (E^(-(1/2) x (a x +
Sqrt[-4 b + a^2 x^2])) (Sqrt[-4 b +
a^2 x^2] (C[1] - E^(x Sqrt[-4 b + a^2 x^2])
C[2])- a x (C[1] + E^(x Sqrt[-4 b + a^2 x^2]) C[2])))
$$f(x)=\frac{e^{-\frac{1}{2} x \left(\sqrt{a^2 x^2-4 b}+a x\right)} \left(\sqrt{a^2 x^2-4 b}
\left(c_1-c_2 e^{x \sqrt{a^2 x^2-4 b}}\right)-a x \left(c_2 e^{x \sqrt{a^2 x^2-4
b}}+c_1\right)\right)}{b}$$
First order, three variables: $(\partial_x+z \partial_z)f(x,y,z)=0$:
In[151]:= symbolicDSolve[
D[f[x, y, z], {x, 1}] + D[f[x, y, z], {y, 1}],
f[x, y, z]] // FullSimplify
Out[151]=
Inactive[Sum][(
x^n252 y^n253 z^
n254 Pochhammer[1 - n252 - n253, -1 + n252] C[1][n254][
n252 + n253])/
Gamma[1 + n252], {n252, 0, \[Infinity]}, {n253,
0, \[Infinity]}, {n254, 0, \[Infinity]}]
$$f(x,y,z)=\underset{\text{n252}=0}{\overset{\infty }{\sum}}\underset{\text{n253}=0}{\overset{\infty }{\sum}}\underset{\text{n254}=0}{\overset{\infty }{\sum}}\frac{x^{\text{n252}}
y^{\text{n253}} z^{\text{n254}} (-\text{n252}-\text{n253}+1)_{\text{n252}-1}
c_1(\text{n254})(\text{n252}+\text{n253})}{\Gamma (\text{n252}+1)}$$
Among doubly infinite solution space, let us choose the following one: $$c_1(n)(k)=\frac{1}{n!}e^{-k}$$. Then, we get
In[165]:=
Sum[Sum[(x^n249 y^(k - n249) z^n251 Pochhammer[1 - k, -1 + n249])/
Gamma[1 + n249] Exp[-k]/(n251 !), {n251, 0, Infinity}, {k, 0,
Infinity}], {n249, a, Infinity}] // FullSimplify
Out[165]= E^z x^a y^-a Gamma[a] Hypergeometric2F1Regularized[1, a,
1 + a, x/y]
which reads
$$f_a(x,y,z)=e^z x^a y^{-a} \Gamma (a) \, _2\tilde{F}_1\left(1,a;a+1;\frac{x}{y}\right)$$
We summed over one of the parameters not from 0 but from $a$ to $\infty$ as the summation is not convergent at the origin around which we have made our series expansion within the code.
If we check the action of our differential operator on our result, we see that
$$(\partial_x+z \partial_z)f_a(x,y,z)=e^{z}x^{-1+a}y^{-a}\sim f_{a-1}(x,y,z)$$ in terms of its divergence around $x\sim 0$. As it stands, it satisfies the differential equation in leading order in regions away from the origin. We can also improve our code such that series expansion point can be chosen manually to avoid singular points.
Mathematicais only the lineDSolve[f[x, y] + D[D[f[x, y], y], x] == 0, f, {x, y}]. The programm doesn't know $f(x,y) = sin(x,y)$. $\endgroup$