Use finite difference method for second ode

I have the following problem. Please help me to solve it:

y" + x^2 * y = sin(\pi * x) 

the boundry conditions are

y(0) = y'(1) =0

We have 20 equals-sub intervals of length h.

Plot the exact solution and the approximate solution.

After solving, I got:

My code to solve this problem is:

%% Boundary Conditions:
y0=0; 
%% Mesh Discretization
N=20; % number of subintervals  
x0=0; xn=1;  % define the boundary nodes
h = (xn-x0)/N; % define the step size (intervals width)
x=(x0:h:xn)';% define the nodes of the whole interval
x_interior = x(2:end-1); % define the interior nodes
m=size(x,1)-2; % define the size of the system (Note since we have the vlaue
             % we subtract 2 from the full size of x
%% Define the matrix M 
M=sparse(m,m);
M(1,1)=1
for i=2:m-1 
    M(i,i)= -2/(h^2) + x_interior.^2; % difine the main diagonal 
end

for i=2:m-1
    M(i-1,i)= (1/h^2)  ; % define the upper diagonal
    M(i, i-1) = (1/h^2); % define the lower diagonal
end

    
    M(m,m)= 3/(2*h);
    M(m,m-1)= -4/(2*h);
    M(m, m-2)= 1/(2*h);

%% define the load vector F 
F= sin ((x_interior.)* pi);
%% Solve the system 
y=M\F; 
% print the value of the solution 
fprintf (['y1= ', num2str(y(1)),'\n', 'y2= ', num2str(y(2)),'\n', 'y3= ', num2str(y(3))])
% plot the solution 
plot(x, [y0; y; yn], 'r');
xlabel ('x')
ylabel( 'y')
hold on 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%% Use matlab bvp4c function 
solinit=bvpinit(x, @guess);
sol=bvp4c (@bvpfcn, @bcfn , solinit);
ref_sol=sol.y(1,:);
plot(sol.x, ref_sol, 'b--')
lg=legend({'FD solution','Reference solution'},'Location','south','FontSize',12);

Where @guess is

function g=guess (x)
g = [0; 0];
end

and @bvpfcn is

function [dydx]= bvpfcn(x,y)
dydx=[y(2); (-x.^2)*y(1)+sin (x.*pi)];
end

and @bcfn is

function res = bcfn (ya, yb)
s=0;  % the value at the end boundary 
res = [ ya(1); yb(1)-s];
end 

I really appreciate your help.