This notebook checks finite difference for second derivative, with uneven position of the x-values.

The first cell
Put vertices at xM = -1,  x0 = 0, and xP = a (with 0 < a <=  1).  These stand for x-minus, x-zero, and x-plus.
Fit a quadratic to uM, u0, and uP.  
The fitted function is f[x]/.soln.
The estimate for u’’ is 2 c2.

Clear[a, f, x, c0, c1, c2, eqns, uM, u0, uP];
f[x_] := c0 + c1 x + c2 x^2
eqns = {f[-1] ==  uM, f[0] ==  u0, f[a] ==  0};
soln = Part[Solve[eqns, {c0, c1, c2}], 1]  (* There is only one soluton, so take the 1st and only one. *)
f[x]/.soln

The fitted function for any h is f[x/h] /. soln /. uVals

h=1.5;  (* The value of h is used only in the graphs *)
a =.5; (* 0 < a ≤ 1, gives position of the boundary point as a fraction of grid spacing *)
u[x_] := -Sin[(x-a h)];
u[x_]:= Exp[a h] - Exp[x];
(* Only change things above this line.  u[x] can be defined in multiple lines, only the last one counts *)

Clear[f, x, c0, c1, c2, eqns, uM, u0];  (* like above, but don't clear a *)
a = N[a]; (* avoids weirdness for a = 1, the integer *)
f[x_] := c0 + c1 x + c2 x^2
eqns = {f[-1] ==  uM, f[0] ==  u0, f[a] ==  0};
soln = Part[Solve[eqns, {c0, c1, c2}], 1];  (* There is only one soluton, so take the 1st and only one. *)
uVals = {uM -> u[-h], u0 ->u[0], uP -> u[a h]};
fit[x_] = f[x/h]/.soln/.uVals;
linFit[x_] = (c0 + c1 x/h)/. soln/.uVals;
p1 = Plot[{u[x], fit[x]}, {x, -h, a h}, PlotLabel ->"u[x] (blue) and fit[x]"];
p2 = Plot[u[x]- linFit[x], {x, -h, a h}, PlotLabel ->"u[x]- linFit[x]"];
p3 = Plot[fit[x]- linFit[x], {x, -h, a h}, PlotLabel ->"fit[x]- linFit[x]"];
p4 = Plot[u[x]- fit[x], {x, -h, a h}, PlotLabel ->"u[x]- fit[x]"];
(*GraphicsGrid[{{p1, p2}, {p3, p4}}, ImageSize ->600] *)(* I only show 2 of these plots *)
GraphicsRow[{p1, p4}, ImageSize -> 700]
hExponentMax = 15;
Print["Whole cell has a = ", a, ", figures have h = ", h];
Print["u[x] = ", u[x], " for this a and h. In table below, function is shifted to make u[a h] = 0."];
eval = -u''[0]/u[0];
Print["Local eval, -u''[0]/u[0] = ", eval];
data = {};
Do[h= 1/2^hExponent;
u[x_] := -Sin[2(x-a h)];
upp = u''[0];
uVals = {uM -> u[-h], u0 -> u[0]};
evalEst = -2c2/h^2/c0/.soln/.uVals;
uppEst = 2 c2/h^2 /.soln/.uVals;
AppendTo[data, {1/h, evalEst, evalEst - eval, uppEst, (uppEst-upp), (uppEst-upp)/h, (uppEst-upp)/h^3}]
,{hExponent, 1, hExponentMax}]
TableForm[data, TableHeadings->{None, {"1/h", "Local Eval est", "LE error", "u''[0]est", "u''error", "u''error/h", "u''error/h^3"}}]

Results for 0-Dirichlet BCs:
The Local Eigenvalue error is O(1).  I assume this is because the denominator, u(0), is going to 0.  But I don’t really understand this and how the local eigenvalue error relates to the eigenvalue error.

The error of u’’(0) is O(h) for 0 < a < 1 and O(h^3) for a = 1.
The prefactor depends on a.

Now assume that the value of u’ at xP = a h is 0.  We are thinking of xP as a boundary point.
We fit two functions u[x], which can be chosen by the uFlag. Both functions have u’(ah) = 0.

As before, we first do a fit with h = 1, and then get a fit for other h by scaling.

Clear[a, f, x, c0, c1, c2, eqns, uM, u0, uP];
f[x_] := c0 + c1 x + c2 x^2
eqns = {f[-1] ==  uM, f[0] ==  u0, f'[a] ==  0};
soln = Part[Solve[eqns, {c0, c1, c2}], 1]

h= 1; (* The value of h is used only in the graphs *)
a = .5;
uFlag =1;  (* uFlag = 1, 2, or 3 for different functions *)
(* Only change things above this line *)

Which[uFlag == 1, 
u[x_] = (x - a h + 1) Exp[-(x - a h + 1)] (* u''' [ah] ≠ 0 *), 
uFlag == 2,
u[x_] = Exp[-(x-a h)^2] (* u'''[ah] = 0  *),
uFlag == 3,
u[x_] = Cos[x-a h]
]; 
(*
Re-do the quadratic fit with the value of a used in this cell 
*)
Clear[f, x, c0, c1, c2, eqns, uM, u0, uP];
f[x_] := c0 + c1 x + c2 x^2
eqns = {f[-1] ==  uM, f[0] ==  u0, f'[a] ==  0};
soln = Part[Solve[eqns, {c0, c1, c2}], 1];
ub = u[a h];
uppVal= u''[0];
uVals = {uM -> u[-h], u0 -> u[0]};
fit[x_] = f[x/h]/.soln/.uVals;
linFit[x_] = (c0 + c1/h x)/. soln/.uVals;
p1 = Plot[{u[x], fit[x]}, {x, -h, a h}, PlotLabel ->"u[x] (blue) and fit[x]"];
p2 = Plot[u[x]- linFit[x], {x, -h, a h}, PlotLabel ->"u[x]- linFit[x]"];
p3 = Plot[fit[x]- linFit[x], {x, -h, a h}, PlotLabel ->"fit[x]- linFit[x]"];
p4 = Plot[u[x]- fit[x], {x, -h, a h}, PlotLabel ->"u[x]- fit[x]"];
(*GraphicsGrid[{{p1, p2}, {p3, p4}}, ImageSize -> 600] *)
GraphicsRow[{p1, p4}, ImageSize -> 700] 
hExponentMax = 13;
Print["Whole cell has a = ", a, " and function ", uFlag, ".  Figures have h = ", h];
Print["u[x] = ", u[x], " and u'''[a h] = ", u'''[a h], " for this a and h."];
Print["In table below, function is shifted to make u'[a h] = 0."];
data = {};
Do[h= 1/2^hExponent;
Which[uFlag == 1, 
u[x_] = (x - a h + 1) Exp[-(x - a h + 1)] (* u''' [ah] ≠ 0 *), 
uFlag == 2,
u[x_] = Exp[-(x-a h)^2] (* u'''[ah] = 0  *),
uFlag == 3,
u[x_] = Cos[x-a h]
]; 
ub = u[a h];
uppVal= u''[0];
eval = -u''[0]/u[0];
uVals = {uM -> u[-h], u0 ->u[0]};
fppEst = (2c2/h^2)/.soln/.uVals;
evalEst = -fppEst/u[0];
ubEst = f[a h]/.soln/.uVals;
AppendTo[data, {1/h,evalEst, evalEst - eval, (evalEst - eval)/h, (evalEst - eval)/h^2, ubEst-ub, (ubEst- ub)/h^2}]
,{hExponent, 1, hExponentMax}]
TableForm[data, TableHeadings->{None, {"1/h", "Local Eval est", "LE error", "LE error/h", "LE error/h^2", "ubEst-ub", "(ubEst-ub)/h^2"}}]

Here are the observed results:  In the table 
Local Eval est = u’’[0]/u[0]est- u’’[0]/u[0], and ub = u(ah) is the value of u at the boundary.  Note that the Local Eval depends on h for all functions except function 3 (the eigenfunction cos).

For function 1, where u''' (ah) != 0, the error of LE estimate is O(h), and the error of ub estimate is O(h^2).

For functions 2 and 3, where u’’’ (ah) = 0, the error of LE estimate is O(h^2), and the error of ub estimate is O(h^2).