## A PDE with non-smooth inhomogeneity
Let $\mathcal{L}$ be a second-order, linear, elliptic differential operator acting on $\mathcal{C}^2([0,2]^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.

Put differently, I solve two PDEs which need to be connected at $x=1$:
\begin{align*}
\begin{cases}
\mathcal{L}u(x,y)=0 & \text{for }(x,y)\in[0,1]\times[0,2], \\
\mathcal{L}u(x,y)+x-1=0& \text{for }(x,y)\in(1,2]\times[0,2].
\end{cases}
\end{align*}


Approximating  all partial derivatives by _central_ differences, I get the nine-point stencil
\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*}
Thus, $u$ is the solution to a system of linear equations.

## 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 the **derivative is not smooth** at $x=1$. The above FD scheme works fine with value-matching (the solution $u$ is perfectly continuous) but struggles with smooth-pasting at $x=1$ (the derivative is not smooth).

[![enter image description here][1]][1]

## Question

> ***How do I ensure smooth-pasting with a finite difference scheme at $x=1$?***


----------


Some of my failed attempts include

 - Impose that forward and backward differences at $x=1$ equal each other (= 2nd order central difference is zero).
 - Use higher order approximations around $x=1$ such as $u_{xx}\approx \frac{-u(-2h)+16u(-h)-30u(0)+16u(h)-u(2h)}{12h^2}$ and $u_{x}\approx \frac{u(-2h)-8u(-h)+8u(h)-u(2h)}{12h}$ and central differences everywhere else.
 - Approximating partial derivatives using points only from one side of $x=1$ (i.e. only using either $u(0), u(h), u(2h)$ or instead $u(0),u(-h),u(-2h)$).
 - Imposing that $u_x\approx\frac{u_{i+1,j}-u_{i-1,j}}{2h}$ equals an average of partial derivatives at $1+h$ and $1-h$.

**Note:** This problem arises as part of a larger system of free boundary problems. Thus, it's necessary to solve the PDE numerically. This question is also posted [here](https://math.stackexchange.com/questions/4011384).


  [1]: https://i.sstatic.net/CTTIF.png