Solving cauchy hyperbolic second order pde I'm currently taking a course in partial differential equations. I'm trying to solve the following problem (which is, as far as I can tell, a bit above the level of the course):

$$\begin{align} u_{xy}&=u \\ \Omega &= \{ x \le y , -\infty \lt x \lt \infty \} \\
  u(x,x) &=  \begin{cases} 0   & x<0 \\ x^2 & x \ge 0
 \end{cases} \\ u_y(x,x) &= 0 \\ \end{align} $$
(using subscript notation to denote partial derivatives)
$u(x,y)$ is continuously differentiable in the solution domain
($\Omega$), and analytic outside some curves.

I tried solving by modifying some methods (e.g. d'Alembert's formula) but it didn't work, I could also identify the equation as 1-dimensional Klein–Gordon Equation and found some particular solutions to the equation, but not to my problem.
How can I solve such a problem?
I prefer a general approach for solving similar problems, but a specific solution will be good enough.
 A: From what I can gather from  my old Graduate PDE course, what you want to do is use something called a Riemann function for the PDE.  I refer to my text on the subject, Guenther & Lee, Partial Differential Equations of Mathematical Physics and Integral Equations, Prentice-Hall (1988), Sec. 4-6, pp. 114-121.
A Riemann function is akin to a Green's function for second-order ODE's and some PDE's.  It is a mapping $R:\mathbb{R}^4 \to \mathbb{R}$, and is used as follows.
Consider the following integral:
$$\int_x^y dx' \, \int_{x'}^y dy' R(x,y,x',y') [u_{x'y'}(x',y')-u(x',y')] = 0$$
The integral is over a triangular region above the line $y'=x'$ and thus in the domain of the given PDE.  Through a bit of trickery, we may rewrite the above integral equation as follows.
$$\int_x^y dx' \, \int_{x'}^y dy' [(R u_{x'})_{y'} - (R_{y'} u)_{x'} + (R_{x'y'}-R) u] = 0$$
We may try to solve this integral equation in many ways, but one way is to define $R$ so as to eliminate unknown quantities.  For example, as we do not know $u$ in the integration region, we may require that
$$R_{x'y'}(x,y,x',y') - R(x,y,x',y')=0 \quad x \lt x',  y \lt y'$$
We may also carry out the other integrations and form other requirements on $R$.  I will summarize.  When the following conditions are met:
$$R_{x'}(x,y,x',y) = 0 \quad x \lt x'$$
$$R_{y'}(x,y,x,y') = 0 \quad y \gt y'$$
$$R(x,y,x,y)=1$$
Then
$$u(x,y) = R(x,y,y,y) h(y) - \int_x^y dx' \, R(x,y,x',x') \phi(x') - \int_x^y dy' R_{y'}(x,y,y',y') h(y') $$
where $h(y)=u(y,y)$, $\psi(y) = u_y(y,y)$, and $\phi(y) = h'(x)-\psi(x)$.  The equation defining $\phi$ is a result of a consistency condition: 
$$h'(x) = \phi(x)+\psi(x)$$
Recall we are given $h(x) = x^2 \theta(x)$ and $\psi(x)=0$, $\theta$ being the Heaviside step function.  NB $h'(x) = 2 x \theta(x) + x^2 \delta(x)$.
I refer you to the cited text as to determining $R$ given the above conditions.  The result is
$$R(x,y,x',y') = J_0{\left ( 2 \sqrt{(x'-x)(y-y')} \right )} $$
$$R_{y'}(x,y,y',y')= -\frac{J_1{\left ( 2 \sqrt{(y'-x)(y-y')} \right )}}{\sqrt{(y'-x)(y-y')}} (x+y-2 y')$$
where $J_0$ and $J_1$ are the Bessel functions of the first kind, zeroth and first order respectively.  The solution is then
$$\begin{align}u(x,y)  &= x^2 \,\theta(x) + \int_x^y dy' \frac{J_1{\left ( 2 \sqrt{(y'-x)(y-y')} \right )}}{\sqrt{(y'-x)(y-y')}} (x+y-2 y') y'^2 \, \theta(y')\end{align}$$
ADDENDUM
Here's a plot of $u(x,x+k)$ for $k \in \{0,1,2,\ldots,9\}$:

