calculation of inverse matrix

I need to calculate the inverse matrix of P2 but I dont know my matrix is sparse or not. also I knew that this problem was solved (I saw the results in a book)

but when I raise 'J', (2^(J+1) is the dimention of my space) I get the singular matrix error even for J=3 however the book has solved for J=6, I dont know exactly in this situation It whould be better to use matlab built-in function or there is more effective algorithm for invers matrix calculation.

This code indeed solve a PDE by wavelet method .

Thanks in advance

% non homogeneous Haar wavelet 

close all
clear all
clc

%% wavelet Twaveoluton
co=1/24;
J=4;
M=pow2(J);
M2=2*M;
A=0;
B=1;
dx=(B-A)/M2;
x(M)=A+M*dx;
ksi1(1)=A; ksi2(1)=B; ksi3(1)=B;
ksi1(2)=A; ksi2(2)=x(M); ksi3(2)=B;
H=zeros(4);
P1=H;
P2=H;
P3=H;
P4=H;
for j=1:J
    
    m=pow2(j);
    mu=round(M/m);
    x(mu)=A+mu*dx;
    x(2*mu)=A+(2*mu)*dx;
    ksi1(m+1)=A;
    ksi2(m+1)=x(mu);
    ksi3(m+1)=x(2*mu);
    for k=1:m-1
        i=m+k+1;
        x(2*k*mu)=A+2*k*mu*dx;
        x((2*k+1)*mu)=A+(2*k+1)*mu*dx;
        x(2*(k+1)*mu)=A+(2*(k+1)*mu)*dx;
        ksi1(i)=x(2*k*mu);
        ksi2(i)=x((2*k+1)*mu);
        ksi3(i)=x(2*(k+1)*mu);
    end
end
xc(1)=0.5*(x(1)+A);
for l=2:M2
    xc(l)=0.5*(x(l)+x(l-1));
    c(l)=(ksi2(l)-ksi1(l))/(ksi3(l)-ksi2(l));
end
for i=1:M2
    K(i)=0.5*(ksi2(i)-ksi1(i))*(ksi3(i)-ksi1(i));
    for l=1:M2
        X=xc(l);
        if X<ksi1(i)
            H(i,l)=0;
            P1(i,l)=0; P2(i,l)=0;P3(i,l)=0;P4(i,l)=0;P6(i,l)=0;
        elseif X<ksi2(i)
            H(i,l)=1;
            P1(i,l)=X-ksi1(i);
            P2(i,l)=0.5*(X-ksi1(i)).^2;
            P3(i,l)=(X-ksi1(i)).^3/6;
            P4(i,l)=co*(X-ksi1(i)).^4;
            P6(i,l)=1/(factorial(6))*(X-ksi1(i))^6;
        elseif X<ksi3(i)
            H(i,l)=-c(i);
            P1(i,l)=c(i)*(ksi3(i)-X);
            P2(i,l)=K(i)-0.5*c(i)*(ksi3(i)-X).^2;
            P3(i,l)=K(i)*(X-ksi2(i))+(ksi3(i)-X).^3/6;
            P4(i,l)=co*((X-ksi1(i)).^4-2*(X-ksi2(i)).^4);
            P6(i,1)=1/(factorial(6))*((X-ksi1(i))^6-(1+c(i))*(X-ksi2(i))^6);
        elseif X>=ksi3(i)
            H(i,l)=0;
            P1(i,l)=0; P2(i,l)=K(i);
            P3(i,l)=K(i)*(X-ksi2(i));
            P4(i,l)=co*((X-ksi1(i)).^4-2*(X-ksi2(i)).^4+(X-ksi3(i)).^4);
            P6(i,1)=1/(factorial(6))*((X-ksi1(i))^6-(1+c(i))*(X-ksi2(i))^6+c(i)*(X-ksi3(i))^6);
            
        else
        end
        F(l)=X-1-0.15*X.^2-X.^5./factorial(5)+sin(3/2.*X).*sin(X./2);
        
    end
    
end
%% calculation the solution of the sine Gordon PDE by wavelet expanstion : 1/L^2V"-V''=sin v

L=10;
dt=1e-3;
tfin=30;
tin=0;
S=(tfin-tin)/dt;
%t=zeros(1,S);
L=20;
beta=0.025;
alpha=L/(sqrt(1-L^2*beta^2));
A=zeros(S+1,2*M);
Vzegx=A;
V=A;
Vdotzegx=A;
Vdot=A;
t=zeros(S+1,1);
B=A;
for d=1:S+2
    t(d)=tin+(d-1)*dt;
end
si=(-4*alpha*beta*exp(-alpha*beta.*t))./(1+exp(-2*alpha*beta.*t));
sidot =(4*alpha^2*beta^2*exp(-alpha*beta.*t).*((1-exp(-2*alpha*beta.*t))))./...
    ((1+exp(-2*alpha*beta.*t)).^2);
si2dot=(4*alpha^3*beta^3*exp(-alpha*beta.*t).*(-1+5*exp(-4*alpha*beta.*t+...
    5*exp(-2*alpha*beta.*t)-exp(-6*alpha*beta.*t))))./((1+exp(-2*alpha*beta.*t)).^4);

fi=4*atan(exp(alpha.*t));
fidot=(-4*alpha*beta*exp(-alpha*beta.*t))./(1+exp(-2*alpha*beta.*t));
fi2dot=(4*alpha^2*beta^2*exp(-alpha*beta.*t).*(1+exp(-2*alpha*beta.*t)))...
    ./((1+exp(-2*alpha*beta.*t)).^2);

for c=1:S+1
    %for f=1:M2
    
    P2inv=eye(2*M)/P2;
    if c==1
        V(1,:)=4*atan(exp(alpha.*xc));
        Vzegx(1,:)=(4*alpha^2.*exp(alpha.*xc).*(1-exp(2*alpha.*xc)))/...
            (1+exp(2*alpha.*xc)).^2;
        
        B(1,:)=(1/L^2.*Vzegx(1,:)-sin(V(1,:))-fi2dot(2).*ones(1,2*M)-...
            si2dot(2).*xc);
        BB=B';
        A(1,:)= B(1,:)/P2;
        Vdot(1,:)=0;
        V2dot(1,:)=0;
        Vdotzegx(1,:)=0;
    elseif c>=2
        V(c,:)=0.5*dt^2*A(c-1,:)*P2+V(c-1,:)+dt*Vdot(c-1,:)+...
            (fi(c)-fi(c-1)-dt*fidot(c-1)+(si(c)-si(c-1)-...
            dt*sidot(c-1)).*xc);
        Vdotzegx(c,:)=dt*A(c-1,:)*H+Vdotzegx(c-1,:);
        Vzegx(c,:)=0.5*dt^2*A(c-1,:)*H+Vzegx(c-1,:)+dt*Vdotzegx(c-1,:);
        B(c,:)=(1/L^2*Vzegx(c,:)-sin(V(c,:))-fi2dot(c+1).*ones(1,2*M)-...
            - si2dot(c+1).*xc);
        A(c,:)=B(c,:)/P2;
        Vdot(c,:)=dt*A(c-1,:)*P2+Vdot(c-1,:)+(fidot(c)-fidot(c-1)).*ones(1,2*M)+(sidot(c)-sidot(c-1)).*xc;
    end
end
% dd=[1,100]
% for d=1:length(dd)
%
%  hold on
plot(xc,V(1000,:))
% end
% error=(B(1000,:)/P2)*P2-B(1000,:);
% cond(A(1000,:))