Impose PDE itself as Boundary Condition?
The underlying problem is that the boundary condition you have provided does not assist in evaluating the solution.
Consider that, for your PDE, if $u=f(x,y)$ is a solution, then so is $u=f(x,y)+C$ for some constant, $C$. This constant cannot be determined using your boundary condition.
Indeed, consider that the following expression is a solution for any choice of constants $C_i$: $$ u=C_1+C_2(x-y)+C_3e^{-x}+C_4e^{-y}+C_5e^{-x-y} $$ None of these constants can be determined by the boundary condition as provided.
That being said, there are often implied conditions that can be used to obtain boundary conditions. For example, if you are working in a real-world context, you may need to apply an energy-conservation requirement. Alternatively, boundedness and related concepts often imply some form of boundary condition. For instance, consider $$ xy'+y=0 $$ In the limit as $x\to0$, we get $y=0$ unless $y'$ tends to infinity. Indeed, the general solution is $y=C/x$, and so if $y$ must remain bounded, then $y(0)=0$, and thus $y=0$.
However, this is actually a separate condition on the problem that generates a boundary condition. You cannot simply apply the DE to the boundary, as it provides no additional information, and boundary conditions are applied to provide information that is otherwise missing.
The mistake you make is the claim that you must apply something at the boundary. This is untrue. You must apply a condition that allows you to identify the correct solution (unless you want a family of solutions, that is).
The nature of the condition will depend on the goal of the differential equation. As mentioned above, sometimes you need to apply some energy-conservation condition. Other conditions are possible - minimisation conditions are common. For example, you might want the solution to your problem that has the least amount of variation (that is, minimise $\iint u_{x}^2+u_{y}^2\ dA$).
There is one other category - cases where your choice won't matter. This is much like how, when performing integration by parts, you don't need to include the constant of integration for the "$\int dv$" - it will cancel out anyway. In this situation, pick the simplest boundary condition for the problem - for instance, if $u(0,0)=1$ and $u(0,\infty)=0$, then you might pick $u(0,y)=e^{-y}$.
EDIT: Here, I'll show what your forward-difference application of the PDE at the boundary is actually assuming for the boundary. Assuming you use central difference for first derivatives when not at the boundary, you have... $$ \begin{align*} \frac{u_{1,j}-u_{0,j}}{\Delta x}+\frac{u_{0,j+1}-u_{0,j-1}}{2\Delta y}+\frac{u_{2,j}-2u_{1,j}+u_{0,j}}{(\Delta x)^2}+\frac{u_{0,j+1}-2u_{0,j}+u_{0,j-1}}{(\Delta y)^2}&=0\\ \frac{u_{2,j}-u_{0,j}}{2\Delta x}+\frac{u_{1,j+1}-u_{1,j-1}}{2\Delta y}+\frac{u_{2,j}-2u_{1,j}+u_{0,j}}{(\Delta x)^2}+\frac{u_{1,j+1}-2u_{1,j}+u_{1,j-1}}{(\Delta y)^2}&=0 \end{align*} $$ With a bit of algebra, we can cancel the $u_{2,j}$ terms to get $$\begin{align*} \left(\frac12+\frac{2+\Delta x}{(\Delta y)^2}\right)u_{0,j} - \left(\frac12+\frac2{(\Delta y)^2}\right)u_{1,j} & \\ + \left(\frac{\Delta x}{4\Delta y} - \frac{\Delta x}{2(\Delta y)^2} + \frac{1}{2\Delta y} - \frac{1}{(\Delta y)^2}\right)u_{0,j-1} - \left(\frac{1}{2\Delta y} - \frac{1}{(\Delta y)^2}\right)u_{1,j-1} & \\ - \left(\frac{\Delta x}{4\Delta y} + \frac{\Delta x}{2(\Delta y)^2} + \frac{1}{2\Delta y} + \frac{1}{(\Delta y)^2}\right)u_{0,j+1} + \left(\frac{1}{2\Delta y} + \frac{1}{(\Delta y)^2}\right)u_{1,j+1} & =0 \end{align*}$$ By Taylor expansion around the $(0,j)$ point, letting $u$ be the function at that point and stopping at $O(\Delta)$ (treating both $\Delta x$ and $\Delta y$ as proportional to $\Delta$), we get $$\begin{align*} \frac{\Delta x}{(\Delta y)^2} u - \frac{\Delta x}2 (u_x+u_y+2u_{xy}) - \frac{\Delta x^2}4 u_{xx} + O(\Delta) & =0 \end{align*}$$ where I have left in the terms that come out explicitly... but note that this is equivalent to $$ \frac{\Delta x}{(\Delta y)^2} u + O(\Delta) = 0 $$ And so, if $\Delta x$ and $\Delta y$ are both scaled down in proportion with each other, your boundary condition is effectively $u(0,y)=0$.
However, be careful, as this analysis applies for keeping $\Delta x$ and $\Delta y$ in proportion. I suspect that some issues may arise if you reduce $\Delta x$ without reducing $\Delta y$.