Where is the error in my below code for approximate solution of poisson boundary value problem?

Question

I have to solve the following boundary value problem which is

also it is defined in my Matlab code below, but my code doesn't work. I mean I didn't get the approximate solution of my system.

I want to know where is the problem in my code or just the version of matlab that I have can't compile the kind of function I have used , Thanks

Explanation of method I have used : I have used the finite element method or what we called Galerkin Method based on investigation about assembly matrix and stiffness matrix. I have multiplied the system by weight function which satisfies the boundary condition then I have integrated over elements (integration of elementary matrix over the range ]-1,1[). I have four elementary matrix. For more information about that Method I used please check this paper(page:6,7,8)

Note The error I have got upon the compilation of my code is

The current use of "MatElt2Nd" is inconsistent with it previous use or definition in line 7

Code

function [U] = EquaDiff2(n)
% ----------------------------------
% -d²u/dx² + 6*u   = (-4*x^2-6)exp(x^2)
% u(-1) = 0 u(1)= 0
%----------------------------------
function [Ke, Fe] = MatElt2Nd(x1,x2) 
  % déclaration de la fonction, 
  % function of computing matrix and elementary matrix (assembly matrix)
  % ----------------------------------
  x = [-1:2/n:1]'; % modification d1 of bound d’intégration
  K = zeros(n+1) ;
  F = zeros(n+1,1) ;
  for i = 1:n
    j = i+1;
    t = [i j];
    x1 = x(i);
    x2 = x(j);

    [Ke,Fe] = MatElt2Nd(x1,x2);
    K(t,t) = K(t,t) + Ke;
    F(t) = F(t) + Fe;
  end;
  K(1,:) = [];
  K(:,1) = [];
  F(1) = [];
  U = K\F;
  U = [0.0;U];
  t = 0:0.01:1;
return

%-------------------------------------------
% calculation of matrix  Ke and vector Fe
%-------------------------------------------
function [Ke,Fe] = MatElt2Nd0(x1,x2) 
  % NEWly named nested function is introduced

  Ke1 = 1/(x2-x1)*[ 1 -1     % no modification done
                    -1 1 ] ; % essentiellement que les matrices
  Ke2 =(x2-x1)* [ 2 1        % élémentaires
                 1 2 ] ;
  N = [(x-x2)/(x1-x2) (x-x1)/(x2-x1)] % function of form
  Fe =simple( int(N' * (-4*x^2-6)*exp(x^2) , x, x1, x2) ) % vecteur Fe ;
  Ke = Ke1 + 6*Ke2 ;

return

Edit I have got a general code for that but I can't do changes in the general code to solve my system , Any help ?

General Code

% au'(x)+bu"(x)=0 for 0<=x<=d
% BC: u(0)=0 and u(d)=h
%==============================================================
% ======Example======
% Finding an approximate solution to the following BVP using 4 elements of
% equal length.
% u'(x)-u"(x)=0 : 0<=x<=1
% BC: u(0)=0 and u(1)=1
% Solution:
% >> Galerkin(4,1,-1,1,1)
% ==============================================================

% The output of this program is 
% 1- The approximate solution (plotted in blue)
% 2- The exact solution (plotted in red)
% 3- The percentage error (plotted in magenta)

%=======================Program Begin==========================
function Galerkin(ne1,a,b,d,h) % Declare function
clc % Clear workspace

% Define the Coefficients of the  exact solution
% The Exact solution is : u(x)=C1+C2*exp(-ax/b)
% where C2=h/(exp(-a*d/b)-1)and C1=-C2
C2=h/((exp(-a*d/b))-1);
C1=-C2;

% Define element length
le = d/ne1; 

% Define x matrix
x = zeros (ne1+1,1); % 
for i=2:ne1 +1
  x(i,1) = x(i-1,1)+le;
end

% K1 matrix corresponding to the diffusion term (u"(x))
K1 = (b/le) * [1,-1;-1,1]
% K2 matrix corresponding to the convection term (u'(x))
K2 =  a*[-1/2 1/2;-1/2 1/2]

% Element stiffness Matrix
Ke = K1+K2

% Global stiffness matrix 
%********************Begin Assembly***************************
k = zeros(ne1+1);
for i=1:ne1+1
  for j=1:ne1 +1
      if (i==j)
          if(i==1)
              k(i,j)=Ke(1,1);
          elseif(i==ne1+1)
              k(i,j)=Ke(2,2);
          else
              k(i,j)=Ke(1,1)+Ke(2,2);
          end
      elseif(i==j+1)
          k(i,j)=Ke(1,2);
      elseif(j==i+1)
          k(i,j)=Ke(2,1);
      else
          k(i,j)=0;
      end
  end
end
%********************End Assembly*****************************


%The Global f Matrix
f = zeros(ne1+1,1);
%BC apply u(0) = 0
f(1,1) = 0;
%BC apply u(d) = h
f(ne1+1,1) = h;
% Display the Global stifness matrix before striking row
K_Global=k
%Striking first row (u1=0)
k(1,1) = 1;
for i=2:ne1+1
  k(1,i) = 0;
  k(ne1+1,i) = 0;
end
k(ne1+1,ne1+1) = 1;
% Display the solvable stifness matrix 
K_strike=k
%solving the result and finding the displacement matrix, {u}
u=inv(k)*f
hold on

% ======Calculating Approximate Solution and plotting============
syms X
U_sym=sym(zeros(ne1,1));
dU_sym=sym(zeros(ne1,1));

for i=1:ne1
    N1x=1-((X-x(i))/le);
    N2x=(X-x(i))/le;
    U_X=(u(i)*N1x)+(u(i+1)*N2x);
    U_sym(i)=U_X;  
    dU_sym(i)=diff(U_sym(i));
    subplot(3,1,1)
    hold on
    ezplot(U_sym(i),[x(i) x(i+1)])
    subplot(3,1,2)
    hold on
    % du/dx approximate
    ezplot(dU_sym(i),[x(i) x(i+1)])
    
end