Let us discuss here. Firstly, I would like point out that in your formula of $x_i$, the range of subscript should be $i=1\dots N+1$. In addition, RHS is means the right hand side of your numerical scheme.
For your definition of partitions, the partition points is given by the formula of $x_i$. Then the step is given by $h_i=x_{i+1}-x_{i}$ for $i=1\dots N$ (for the sake of simplicity, I use $h_i$ instead of delta_x). By doing this, the matrix A should be reformed in $h_i$, because the original program you wrote is in the case of uniform mesh. What you should do is to adjust the coefficient delta_x to $h_i$. But you must be care of in this case, the step length is depend on step number. What is more, you also should adjust the RHS of your numerical scheme because in your original code, the RHS, say $f(x_i)$, is given in a uniform partition. But in your singular mesh, the RHS should given in the $x_i$ you defined, which is different from uniform mesh.
In summary, just make two vectors for partition points $x_i$ and step length $h_i$. And then modify the code as I told you, I think you can solve the problem.
A sample code:
clear;
N=10;
ax=0;
bx=1;
px=1-cos(pi/2*(0:N)/N);
h=px(2:end)-px(1:end-1);
exau=@(x)x.^2;
RHSf=@(x)-2+zeros(size(x)); %-u''=f
u=zeros(N+1,1);
u(1)=exau(ax);
u(end)=exau(bx);
A=diag(1./h(2:end)+1./h(1:end-1))-diag(1./h(2:end-1),1)-diag(1./h(2:end-1),-1);
RHS=RHSf(px(2:end-1)').*(h(1:end-1)+h(2:end))'/2;
RHS(1)=RHS(1)+u(1)/h(1);
RHS(end)=RHS(end)+u(end)/h(end);
u(2:end-1)=A\RHS;
plot(px,u-exau(px'));
% hold on;
% plot(linspace(ax,bx,N+1),exau(linspace(ax,bx,N+1)'),'-r');
max(abs(u-exau(px')))
The numerical scheme of 2nd derivative:
Let $x_i$ is the partition point and $h_i=x_{i+1}-x_i$ is the step length. What is more, $u_i$ denote $u(x_i)$ and $u_{i+1/2}$ denote value of $u(x)$ at the middle point of $[x_i,x_{i+1}]$.
\begin{equation}
u_i''=\frac{u'_{i+1/2}-u'_{i-1/2}}{h_{i+1}/2+h_i/2}
\\=\frac{\frac{u_{i+1}-u_i}{h_{i+1}}-\frac{u_i-u_{i-1}}{h_i}}{h_{i+1}/2+h_i/2}
\end{equation}
Just render it in a neat form, you will get the numerical scheme.
Analytic solution :
$$y(x)=-\int \frac{-ax+b+k}{-x^2+x+m} = \frac{a-2b-2k}{\sqrt{4m+1} }\tanh^{-1}\left( \frac{2x-1}{\sqrt{4m+1}} \right) - \frac{a}{2}\ln|-x^2+x+m | +C$$
With condition $y(0)=-y(1)$ :
$y(0)=\frac{a-2b-2k}{\sqrt{4m+1} }\tanh^{-1}\left( \frac{-1}{\sqrt{4m+1}} \right) - \frac{a}{2}\ln|m | +C$
$C=y(0)-\frac{a-2b-2k}{\sqrt{4m+1} }\tanh^{-1}\left( \frac{-1}{\sqrt{4m+1}} \right) + \frac{a}{2}\ln|m|$
$$y(x)=y(0)-\frac{a-2b-2k}{\sqrt{4m+1} }\tanh^{-1}\left( \frac{-1}{\sqrt{4m+1}} \right) + \frac{a}{2}\ln|m|+ +{\frac{a-2b-2k}{\sqrt{4m+1} }\tanh^{-1}\left( \frac{2x-1}{\sqrt{4m+1}} \right) - \frac{a}{2}\ln|-x^2+x+m |} $$
As a consequence :
$y(1)=y(0)-\frac{a-2b-2k}{\sqrt{4m+1} }\tanh^{-1}\left( \frac{-1}{\sqrt{4m+1}} \right) + \frac{a}{2}\ln|m| + {\frac{a-2b-2k}{\sqrt{4m+1} }\tanh^{-1}\left( \frac{1}{\sqrt{4m+1}} \right) - \frac{a}{2}\ln|m|} $
$y(1)=y(0)+2\frac{a-2b-2k}{\sqrt{4m+1} }\tanh^{-1}\left( \frac{1}{\sqrt{4m+1}} \right)$
The second condition $y(1)=-y(0)$ implies :
$-y(0)=y(0)+2\frac{a-2b-2k}{\sqrt{4m+1} }\tanh^{-1}\left( \frac{1}{\sqrt{4m+1}} \right)$
$y(0)=-\frac{a-2b-2k}{\sqrt{4m+1} }\tanh^{-1}\left( \frac{1}{\sqrt{4m+1}} \right)$
$$y(x)=-2\frac{a-2b-2k}{\sqrt{4m+1} }\tanh^{-1}\left( \frac{-1}{\sqrt{4m+1}} \right) + \frac{a}{2}\ln|m|+ +{\frac{a-2b-2k}{\sqrt{4m+1} }\tanh^{-1}\left( \frac{2x-1}{\sqrt{4m+1}} \right) - \frac{a}{2}\ln|-x^2+x+m |} $$
Best Answer
I saw your code there. Just use it, and replace the boundary conditions with the following: $$ u'(0)\approx\frac{u_1-u_0}{\Delta x}, $$ $$ u'(1)\approx\frac{u_M-u_{M-1}}{\Delta x}, $$ where $u_j$ is your numerical solution.
Thus you have: $$ u_0=u_1-u'(0)\Delta x, $$ $$ u_M=u_{M-1}+u'(0)\Delta x. $$
This is a first order approximation. Expand further the Taylor series (using more grid points) to enhance the approximation goodness.
In your specific case: $$ u_0=u_1, $$ $$ u_M=u_{M-1}. $$