two-point boundary value problem with piecewise function coefficients Let $u(x)$ in $(0,1)$ satisfy：$$\frac{d}{dx}(a(x)\frac{du}{dx}) = 1$$ With $u(0) = u(1) = 0$. $a(x)$ is a piece-wise function：$$a(x) = \left\{\begin{aligned}
    1,& & x\in(0,1/2) \\
    2,& & x\in(1/2,1) \\
\end{aligned}\right.$$ It also requires that $u(x)$ and $a(x)\frac{du}{dx}$ be continuous at $x=1/2$.
And the question is to solve this equation using finite-difference methods (use Matlab)
I write the code the textbook way (when $a(x)$ is continuous), but the result is quite strange:

I use the method from this article,P75

Following is function code, please tell me how to fix it if you find something wrong. Thank you!
function u = Ch3_varcof_1d(a,f,N)
%--------------------------------------------------------------------------
% u = Ch3_FD1d(a,b,c,f,N,method)
% 
% compute the numerical solution of the 1D example (h=1/N)
%   d/dx(a(x)du/dx) = f  in (0,1)
%   u(0) = u(1) = 0
% 
%               2021-05-18
%--------------------------------------------------------------------------

h = 1/N; % Step-length
x_list = 0:1/N:1; % Discrete-Points in interval
l1 = ones(N-2,1); % Downside diagonal
l2 = ones(N-1,1); % Middle diagonal
l3 = ones(N-2,1); % Upside diagonal

for i=2:N-1
    l1(i-1) = a((x_list(i)+x_list(i+1))/2);
end

for i=1:N-1
    l2(i) = a((x_list(i+1)+x_list(i+2))/2) + a((x_list(i)+x_list(i+1))/2);
end

for i=2:N-1
    l3(i-1) = a((x_list(i+1)+x_list(i+2))/2);
end

A = (diag(l1,-1)-diag(l2)+diag(l3,1))./(h^2);

if isa(f,'function_handle') % Get the function F vector
    x = (0:h:1)';
    F = f(x(2:end-1));
else
    F = f*ones(N-1,1);
end

u = A \ F;
u = [0; u; 0];

 A: You don't need any specialised method to solve this problem... To make sure that $x=\frac 12$ is one of the nodes, lets use $n = 2^k$ subintervals ($2^k+1$ points). This way the central node has index $m = \frac n2 +1$. Since $a(x)$ is piecewise constant, at every node except from $m$ the equation to be satisfied is $
u''(x) = \frac{1}{a(x)}$, yielding the linear equations
$$
u_1=0
$$
$$
u_{i+1}-2 u_i +u_{i-1} = \frac{h^2}{a(x_i)}, \quad i= 2,\cdots, n-1, \quad i \ne m
$$
$$
u_n=0
$$
Regarding the $m^{th}$ equation, the continuity of $a(x) u'(x)$ yields
$$
\frac{u_m-u_{m-1}}{h} = 2 \times \frac{u_{m+1}-u_m}{h}.
$$
so, you get the standard discretisation of the second derivative, with one modified equation. The system matrix for $n=2^3$ would be
$$
\left(
\begin{array}{ccccccccc}
 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 1 & -2 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & 1 & -2 & 1 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 1 & -2 & 1 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & -1 & 3 & -2 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 1 & -2 & 1 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 1 & -2 & 1 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 1 & -2 & 1 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\
\end{array}
\right)
$$
And graphically, the solution is

If you don't enforce the continuity of $a(x) u'(x)$, you get nicer solutions, which don't present the small "hiccup" at $x=\frac 12$.
