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.
Best Answer
This is called a singular boundary-value problem. Matlab can handle some singular BVPs (look at the documentation for
bvp4cand theSingularTermoption in bvpset) so you need to bring your equation in the form that Matlab can handle.Added later: I've never used this option before, but here is how I would start. When bringing the equation to first-order form, one normally introduces a vector $(F_1,F_2)$ with $F_1(K) = F(K)$ and $F_2(K) = F'(K)$ so that the equation becomes $$ \begin{align} F_1' &= F_2, \\ F_2' &= \frac{3.45}{0.02K^2} F_2 + \frac{0.08}{0.02K^2} F_1 + \frac{3.45}{0.02K^2}. \end{align} $$ As you noticed, Matlab requires the equation to be of the form $y' = \frac1x Sy + f(x,y)$ ...
Added yet later: I initially thought that one could get rid of the $1/K^2$ terms by using $F_1(K) = F(K)$ and $F_2(K) = KF'(K)$ but as Qiang Li noted (in more polite terms) I made a stupid mistake. In fact, in the original equation the boundary point $K=0$ is an irregular singular point, while in the equations that Matlab can handle the boundary point is a regular singular point. My uneducated guess would be that the algorithm implemented in Matlab won't be able to handle your equation. Sorry.