Problem with Singular Matrices

I am working on a crank-nicolson scheme for solving the Allen-Cahn equation. Below is the code that I have written, but I keep on getting the warning of matrix is singular to working precision. I have also attached the paper scheme on which the code is based on for reference.

Many thanks in advance, any help on getting this program running will be greatly appreciated.

Here is my code:

clear all

epsilon=10^(-3);
xL=-5;xR=-xL;
yD=-5;yU=-yD;
L=xR-xL; %the length of domain
tau=0.02;    %time step "k"
h=0.2;      %spatial step
%x=(xL+h):h:(xR-h); y=(yD+h):h:(yU-h);
x=xL:h:xR; y=yD:h:yU;
r=tau/h^2;
N=L/h; T=10;
Utemp=zeros(N+1,N+1);
x0=floor((N+1)/2);y0=x0;
for i=1:N+1
    for j=1:N+1
        if  ((i-x0)^2+(j-y0)^2)<=5^2
            Utemp(i,j)=1;
        end
    end
end
%mesh(Utemp)
U=reshape(Utemp',(N+1)^2,1);
%U=sparse(U);
%mesh(U)
%B1=diag((2+4*r)*ones(N-1,1))+diag(-r*ones(N-2,1),1)+diag(-r*ones(N-2,1),-1);
%B2=diag((2-4*r)*ones(N-1,1))+diag( r*ones(N-2,1),1)+diag( r*ones(N-2,1),-1);

K1 = diag([0,-r*ones(1,N-1),0]);
K11 = diag([0,r*ones(1,N-1),0]);
for t=0:tau:T

      %  for the left parts
      A = kron(diag([ones(1,N-1),0],-1),K1) + kron(diag([0,ones(1,N-1)],1),K1);
      %  for the right parts
      A1 = kron(diag([ones(1,N-1),0],-1),K11) + kron(diag([0,ones(1,N-1)],1),K11);
      for i=2:N
          Ui = Utemp(i,:);
          K2 = diag([1,(2+4*r-tau*epsilon^(-2)*(1-Ui(2:N).^2)),1]) + diag([-r*ones(1,N-1),0],-1) + diag([0,-r*ones(1,N-1)],1);
          A = A + kron(diag([zeros(1,i-1),1,zeros(1,N+1-i)]),K2);

          K22 = diag([1,(2-4*r+tau*epsilon^(-2)*(1-Ui(2:N).^2)),1]) + diag([r*ones(1,N-1),0],-1) + diag([0,r*ones(1,N-1)],1);
          A1 = A1 + kron(diag([zeros(1,i-1),1,zeros(1,N+1-i)]),K22);

      end
       A = kron(diag([1,zeros(1,N-1),1]),eye(N+1));
       A1 = kron(diag([1,zeros(1,N-1),1]),eye(N+1));

       U=(A\A1)*U;
       Utemp=reshape(U,N+1,N+1);
       Utemp=Utemp';
       pause
       mesh(Utemp)
  end
  %U=reshape(Utemp,N+1,N+1);
  mesh(Utemp)

[Walter Roberson]

Attachment did not make it.

[Jamil Villarreal]

Attachment should be available now.

Answer 1

[Walter Roberson]

You have

       A = kron(diag([1,zeros(1,N-1),1]),eye(N+1));
       A1 = kron(diag([1,zeros(1,N-1),1]),eye(N+1));

       U=A\A1*U;

Your A and A1 are the same and are rank 2*(N+1) but size (N+1)^2 x (N+1)^2. You then \ those identical low-rank matrices. That is the cause of the singularity warning.

Remember that A\A1*U is (A\A1)*U not A\(A1*U)

If somehow A and A1 were not singular, then because they are identical, A\A1 would be the numeric approximation of the identity matrix.

[Jamil Villarreal]

Thank you for pointing that out. I did include the parenthesis as suggested, but I still get the warning of singular matrices.

[Walter Roberson]

I am not surprised; if I remember my algebra correctly, matrix multiplication cannot increase the rank.

[Jamil Villarreal]

So you're saying that the size of both A and A1 are causing the singularity, correct? If so, then how could I fix this issue?