clear;
clf;
 
N=50;
L=9.53;
 
% % % %movieeee
% % % set(gca,'nextplot','replacechildren');
% % % aviobj = avifile('222','quality',100)

ppp=0;

for i=1:N+1
    x(i)=(i-1)*L/N;
    y(i)=(i-1)*L/N;
end

T(:,:)=0;
% for i=(1:N);
%     for j=(1:N);
%         if (i-(N/2))^2+(j-(N/2))^2<(400);
%         T(i,j)=100;
%         end
%     end
% end
% for i=(1:N);
%     for j=(1:N);
%         if (i-(N/2))^2+(j-(N/2))^2<(100);
%         T(i,j)=0;
%         end
%     end
% end

%initial conditions
%for i=0:N
    %for j=0:N
        %T(i+1,j+1)= min([i,N-i,j,N-j])*4;
        
    %end
%end

T(1,1:N+1) = 0;
T(N+1,1:N+1) = 0;
T(1:N+1,1) = 0;
T(1:N+1,N+1) = 0;

Tnew=T;

xx=0:L/N:L;
yy=0:L/N:L;
[X,Y] = meshgrid(yy,xx);
surf(X,Y,T);
view([+35 154])
% pause

f=0;

n=.25;
for iter=1:150
    if f<N
        f=f+1;
    end
    
    for i=2:N
        for j=2:N
            Tnew(i,j)= n*(T(i+1,j)+T(i-1,j)+T(i,j-1)+T(i,j+1)-4*T(i,j))+T(i,j);
        end
    end
    for i=20:30
        for j=2:f;
            Tnew(i,j)=Tnew(i,j)+1;
        end
    end
    
    for v=1:N
        Tnew(1,v)=Tnew(2,v);
        Tnew(50,v)=Tnew(49,v);
        Tnew(v,1)=Tnew(v,2);
        Tnew(v,50)=Tnew(v,49);
    end
    
    T=Tnew;
    
    surf(X,Y,T);
    axis([0 L 0 L 0 100]);
    view([+35 154]);
    pause (0.1);
    
% % %     if mod(iter,10) == 1
% % %         aviobj = addframe(aviobj,getframe);
% % %     end
%     saveas(gcf, ['frame',sprintf('%03d',ppp),'.png'])
%     ppp=ppp+1;

end 

Tnew=T;

for iter=1:300
    for i=2:N
        for j=2:N
            Tnew(i,j)= n*(T(i+1,j)+T(i-1,j)+T(i,j-1)+T(i,j+1)-4*T(i,j))+T(i,j);
        end
    end

    for v=1:N
        Tnew(1,v)=Tnew(2,v);
        Tnew(50,v)=Tnew(49,v);
        Tnew(v,1)=Tnew(v,2);
        Tnew(v,50)=Tnew(v,49);
    end
    
    T=Tnew;
    
    surf(X,Y,T);
    axis([0 L 0 L 0 100]);
    view([+35 154]);
    pause (0.1);
    
% % %     if mod(iter,10) == 1
% % %         aviobj = addframe(aviobj,getframe);
% % %     end
%     saveas(gcf, ['frame',sprintf('%03d',ppp),'.png'])
%     ppp=ppp+1;

end

% % % aviobj = close(aviobj);