
clear

n = 51;
X = linspace(-3,3,n);
Y = linspace(-3,3,n);
h = X(2)-X(1);

c = 2;
d = -1;
e = 1.2;
f = 0.75;
k = 1;
shift = 10;

u0 = @(x,y) 1-max(sign((c^2-(x/e-d).^2-(y/f).^2).^2 + k*(c+d-x/e).^3-shift),0);

for i=1:n
    for j=1:n
        Z(i,j) = u0(X(i),Y(j));
    end
end
mesh(Y,X,Z)

T = (-2*diag(ones(n,1))+diag(ones(n-1,1),1)+diag(ones(n-1,1),-1));

A = sparse(1/h^2*(kron(T,eye(n,n))+kron(eye(n,n),T)));

U(:,1) = Z(:);
tau = 0.01;
T = 0:tau:3;

for i=1:length(T)-1
    U(:,i+1)=(eye(n^2,n^2)-tau/10*A)\U(:,i);
end

for i=1:length(T)
    mesh(X,Y,reshape(U(:,i),n,n));
    drawnow
    set(gca,'xlim',X([1 end]),'ylim',Y([1 end]),'zlim',[0 1]);
    M(i) = getframe;
end
movie(M)
