Numerical Code for 2D -Transient heat conduction equation

%% Given Values
clear all;
L=1;                     %edge length
A=0.01;                  %area under consideration
alpha=9.7e-05;           %thermal diffusivity aluminium
n=input('grid size');    %enter grid size
dx=L/n;
dy=dx;
dt=dx/100;
x=linspace(0,1,n);
y=linspace(1,0,n);
time=input('enter time');%enter time
t=zeros(n,n);
t(:,:)=200;              %initial condition

%% Linear boundary condition
t(1,:)=linspace(300,320,n);
t(:,1)=linspace(300,330,n);
t(n,:)=linspace(330,270,n);
t(:,n)=linspace(320,270,n);
%% solving f.d.e
convrg=1;
ts=0;
tol=1e-06;
while ts<=time && convrg>=tol     %convergence criterion
    tpast=t;
    for i=2:n-1 
            for j=2:n-1;
                 t(i,j)=((tpast(i+1,j)-2*tpast(i,j)+tpast(i-1,j))/dx^2 + (tpast(i,j+1)-2*tpast(i,j)+tpast(i,j-1))/dy^2)*alpha*dt+tpast(i,j);
            end
    end
            ts=ts+dt;
convrg=max(max(abs(t-tpast)));
end

%% figures
figure(1);
contour(x,y,t);title('temperature isoline');
xlabel('x position in meter');ylabel('y position in meter');
box on; grid on;
legend('temp in kelvin');colormap('hot');
q=colorbar;ylabel(q,'temperature in kelvin');

figure(2);
pcolor(x,y,t); shading interp;box on; grid on;title('temperature distribution');xlabel('position in meter x-axis');ylabel('postion in meter y-axis');colormap('hot');
h=colorbar;ylabel(h,'temperature in kelvin');

figure(3);
plot(t(1:n,n/2),y);title('tempearture (x=0.5m,y)');xlabel('temperature in K');ylabel('y in mt');legend('x=0.5m');box on; grid on;