N    = 30;
h    = 1/N;
N1   = N+1;
N1N1 = N1*N1;
N4   = round(N/4);

M  = 20;
t  = 1/M;

r  = 1.0;
k  = 1.0;
d  = 1.0;

V    = [zeros(N4,1);ones(N1-2*N4,1);zeros(N4,1)];
V    = kron(V,V');
V    = reshape(V,N1N1,1);

D0 = -eye(N1);
D1 = diag([4*ones(N,1);3])-diag(ones(N,1),+1)-diag(ones(N,1),-1);
D  = kron(eye(N1),D1)+kron(diag(ones(N,1),+1),D0)+kron(diag(ones(N,1),-1),D0);

D2 = diag(ones(N1,1))-diag(ones(N,1),-1);
F  = kron(eye(N1),D2);

B = (d/h^2)*D + (k/h)*F + r*speye(N1N1);

A = speye(N1N1) + t*B;
D = diag(diag(A));
[L,U] = lu(A);

h1 = figure(1); close(h1);
h1 = figure(1); set(h1,'Position',[10 10 300 300]);

method = 1;

for n=1:M
    switch method
        case 1
            if (k == 0)
                V = pcg(A,V,1.0e-6,100000,D);
            else
                V = pcg(A'*A,A'*V,1.0e-6,100000,D);
            end
        case 2
            V = L\V;
            V = U\V;
        case 3
            V = A\V;
    end
    surf(reshape(V,N1,N1)); axis([1 N1 1 N1 0 1]); view([1 1 1]);
    drawnow; pause(0.1);
end
