**Domain for X an Y**: I solve a system of diff equations to get the 1st, 2nd, and 3rd moment of X and 1st and 2nd moment of Y. Then i get the equivalent log-normal displaced process for X and the gaussian one for Y.**Differential operator for X**: I use a tridiagonal operator for X with a first order upwind scheme for the first order derivative (this does not significant affect the results)**Differential Operator for Y**: In the Y-direction i apply a change of variable as described in Piterbarg and Andersen's book. I really solve in a variable where the mean of y is substracted. I have implemented two different discretizations schemes for the derivative in Y.

1. Upwind Scheme: Use of two point approx for the derivative according to the sign of the Y-drift.2. 5 point approximation: i fit a forth order polinomial and derive that polinomial**Boundary Conditions**: I assume that second derivatives (in both directions) go to zero at the boundaries.

- For positive mean reversions (or moderate negative ones) both schemes (1st order upwind scheme and five-point approximation) behaves OK. I have tried up to a 30y-20y Swap.
- For negative mean reversions, results start worsening: For -10% values (for mean reversion) 1st order upwind scheme start getting unstable for IRS starting in 20y time. The five point discretization works much worse for negative values of the mean reversion. With the latter, while killing the volatility (and keeping the XY domain fixed) i recover Values for the IRS (at time 0) extremely big (in absolute values) at both extremes (I am not sure if in this case the solver for the penta-diagonal system migh not be working properly for high values of the npv) .

I would highly appreciate any comment or suggestion that helps me to find where the problem is...