Let $\mathcal{L}$ be a second-order, linear, elliptic differential operator acting on $\mathcal{C}^2([0,\infty)^2)$.

I'm numerically solving the inhomogeneous PDE
\begin{align*}
\mathcal{L}u(x,y)+(x-1)^+=0,
\end{align*}
where $(\cdot)^+$ denotes the positive part.
## Finite Difference Scheme ##
Approximating  all partial derivatives by _central_ differences, I get a FD scheme, $Au=b$, whose solution gives me the grid points $u_{i,j}$. The FD scheme looks like
\begin{align*}
c_1 u_{i-1,j-1} + c_2 u_{i,j-1} + c_3 u_{i+1,j-1} + c_4 u_{i-1,j} + c_5 u_{i,j}&\\
 + c_6 u_{i+1,j} + c_7 u_{i-1,j+1} + c_8 u_{i,j+1} + c_9 u_{i+1,j+1} + (x_i-1)^+ &=0,
\end{align*}
where the $c_k$ are coefficients (independent of $x$ but dependent on $y$). 

## Problem ##
Plotting the solution $u$, it all looks fine and perfect. However, a plot of $\frac{\partial u}{\partial x}$ as a function of $x$ shows that there's a **small kink** at $x=1$. To find an analytical solution in the one-dimensional case, I would impose value-matching and smooth-pasting conditions at $x=1$. The above FD scheme seems to work fine with value-matching (the solution $u$ looks appears perfectly continuous) but struggles with smooth-pasting.

**Question:** How do I ensure smooth-pasting with the FD scheme? I tried to impose that forward and backward differences at $x=1$ equal each other (= 2nd order central difference being zero) but it didn't do the trick.

**Note:** This problem arises as part of a larger system of free boundary problems. Thus, it's necessary to solve the above PDE numerically.