clear
clf
f=0;
N=100;
L=10;
T=zeros(N+1,N+1);
maxX = 101;
maxY = 101;

    x0=maxX/2
    y0=maxY/2
    R=maxX/2.5;
for i=1:maxX;
    for j=1:maxY;
        dx=i-x0+1;
        dy=j-y0+1;
        r=hypot(dx,dy);
        if (r<R);
            T(i,j)=120;
        end
    end

      
end
        x0=maxX/2
        y0=maxY/2
        R=maxX/4;
        for i=1:maxX;
        for j=1:maxY;
            dx=i-x0+1;
            dy=j-y0+1;
            r=hypot(dx,dy);
            if (r<R);
                T(i,j)=260;
            end
        end
        end

x=zeros(1,N+1);
y=zeros(1,N+1);
 
k=111;
dx=3;
dt=.01;
dy=3;
 
for i=1:N+1
    x(i)=dx*(i-1);
    y(i)=dy*(i-1);
end
 
Numit=20000;
 
for q=1:Numit
 
Tnew=T;
 
        for i=2:N

        for j=2:N;
            
            Tnew(i,j)=(k)*(dt/dx^2)*(T(i+1,j)-2*T(i,j)+T(i-1,j) + T(i, j+1) -2*T(i,j) + T(i, j-1))+ T(i,j);
        end

        end

        
 T=Tnew;
    if mod(q+9,80)==0
        f=f+1;
        xx=0:L/N:L;
        yy=0:L/N:L;
        [X,Y] = meshgrid(yy,xx);
        surf(X,Y,T)
        axis([0 L 0 L 0 300]);
        view(180,90)
        %pause
       %saveas(gcf, ['wire_heat_',sprintf('%03d',f),'.png'])
        pause (0.01)
    end
    
    fprintf('%number of iterations have been done/n',q)
  end

