I changed boundary condition as you suggested.

My FDM model and QMC model yield very close result!
s=95,K=90,H=100,t=0.1,sigma=0.2,r=0,monitorfrequency=20
QMC=1.952 (50,000 trials)
FDM=1.951 (2000 steps)

I still have no idea about smin and smax by the way. In this example I just randomly pick smin=50 and smax=150 (so one node in price axis lies exactly on the barrier)

I changed boundary condition as you suggested.

My FDM model and QMC model yield very close result!
s=95,K=90,H=100,t=0.1,sigma=0.2,r=0,monitorfrequency=20
QMC=1.952 (50,000 trials)
FDM=1.951 (2000 steps)

I still have no idea about smin and smax by the way. In this example I just randomly pick smin=50 and smax=150 (so one node in price axis lies exactly on the barrier)

I used the BKG barrier correction formula with the same data as above. I get
m = 20
S, C
====================
94.95917 1.994772
94.97934 1.992696
94.99951 1.990581
95.01967 1.988427
95.03984 1.986233

My FDM model and QMC model yield very close result!
They could both be a bit off?

I changed boundary condition as you suggested.

My FDM model and QMC model yield very close result!
s=95,K=90,H=100,t=0.1,sigma=0.2,r=0,monitorfrequency=20
QMC=1.952 (50,000 trials)
FDM=1.951 (2000 steps)

I still have no idea about smin and smax by the way. In this example I just randomly pick smin=50 and smax=150 (so one node in price axis lies exactly on the barrier)

I used the BKG barrier correction formula with the same data as above. I get
m = 20
S, C
====================
94.95917 1.994772
94.97934 1.992696
94.99951 1.990581
95.01967 1.988427
95.03984 1.986233

My FDM model and QMC model yield very close result!
They could both be a bit off?

Yes they were and so can the BKG correction it seems. The correct price here is 1.9602906. Is it hard to get the MOL working for discrete barriers? And what would its stability characteristics be, would it handle the repeatedly imposed discontinuities as gracefully as the fully implicit setup?

By the way, you haven't confirmed your time scheme, I've been assuming all along it's the (fully) implicit method (thus 1st order in time, thus the above extrapolation in time). I've also been assuming that your spatial (asset) discretization is using central differences so it's second order. And that's why you cannot just double both NS and NT together and do the RE, because you have different orders in asset and time. RE depends on the order.

This is my result.
I pick smin =0.4*H and smax=1.2*H so that there is no saw tooth.
Best I can get is 1.9608 with Payoff-averaging 2,000-step FDM.

That convergence graph looks right. That's already a huge improvement over your "naive" implementation. Now there are probably two reasons why you are still a little far off the exact result. One is that you are using the implicit scheme (are you?) which is first order in time. My result comes from a second order scheme. So if you are using implicit, first thing you should try is Richardson Extrapolation which in this case means calculate first a fine price as cuo(NS, NT) and then a crude one as cuo(NS, NT/2) and find the 2nd order (in time) correct price as 2*fine_price - crude_price. I am assuming you're using central differences otherwise so you're 2nd order in space.
The other reason why you may be a little off is that the payoff-averaging (or intermediate solution smoothing) you used is the simplest one. My code is using something slightly more fancy, but we may get to that later, first try the RE and see what happens.
P.S. I didn't understand your "so that there is no saw tooth." comment. We know there should be because the solution is allowed to become non-zero above the barrier in between monitoring dates, no?

I still see 2000x2000 result more accurate than 2000x40+2000x20 result. I don't know why.

This is my result.
I pick smin =0.4*H and smax=1.2*H so that there is no saw tooth.
Best I can get is 1.9608 with Payoff-averaging 2,000-step FDM.

That convergence graph looks right. That's already a huge improvement over your "naive" implementation. Now there are probably two reasons why you are still a little far off the exact result. One is that you are using the implicit scheme (are you?) which is first order in time. My result comes from a second order scheme. So if you are using implicit, first thing you should try is Richardson Extrapolation which in this case means calculate first a fine price as cuo(NS, NT) and then a crude one as cuo(NS, NT/2) and find the 2nd order (in time) correct price as 2*fine_price - crude_price. I am assuming you're using central differences otherwise so you're 2nd order in space.
The other reason why you may be a little off is that the payoff-averaging (or intermediate solution smoothing) you used is the simplest one. My code is using something slightly more fancy, but we may get to that later, first try the RE and see what happens.
P.S. I didn't understand your "so that there is no saw tooth." comment. We know there should be because the solution is allowed to become non-zero above the barrier in between monitoring dates, no?

Oh sorry, somehow, I missed this reply. So you have repeatedly reply me everything.

First set
I use smin=0 and smax=that formula that I come up myself=125.
Well I just want to make sure that barrier lies on the grid at N=0.8N.
Second set
I use smin=0.8*H=80 and smax=105

First set
I use smin=0 and smax=that formula that I come up myself=125.
Well I just want to make sure that barrier lies on the grid at N=0.8N.
Second set
I use smin=0.8*H=80 and smax=105

A OK.Your Smin is way too high in this case, just use 20 for now and Smax at 120 you still should be able to get very accurate results with 2000 steps.