Plotting 2D Temperature Distribution

Hi, I am having a problem plotting 2d Temperature Distribution on a geometry.

I have solved it using Finite difference method,the temperatures of all the nodes.

That is the proble, i choose the mesh size of 0.05cm, i just need to draw and plot the temperatures obtained from the code the code is as below

clc
clear all

Tinf=295;
Tc=3350;

h_a=5;          % air h
h_c=8800;       % coolant h
h_cc=2350;      % combustion chamber h

l=0.0005;      % mesh size

k=14.9;         % Thermal conductivity 304 Steel

h_a_f=h_a*l/k;
h_c_f=h_c*l/k;
h_cc_f=h_cc*l/k;

Ta=295;
Tc=183;
Tcc=3350;


% The grid has 31 nodes,

t1=[-2-h_a_f,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
t2=[1,-3-h_a_f,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
t3=[0,1,-3-h_a_f,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
t4=[0,0,1,-3-h_a_f,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
t5=[0,0,0,1,-2-h_a_f,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
t6=[1,0,0,0,0,-3,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
t7=[0,1,0,0,0,1,-4,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
t8=[0,0,1,0,0,0,1,-4,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
t9=[0,0,0,1,0,0,0,1,-4,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
t10=[0,0,0,0,1,0,0,0,1,-3,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
t11=[0,0,0,0,0,1,0,0,0,0,-3-h_c_f,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
t12=[0,0,0,0,0,0,1,0,0,0,1,-4-h_c_f,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
t13=[0,0,0,0,0,0,0,1,0,0,0,1,-3-h_c_f,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
t14=[0,0,0,0,0,0,0,0,1,0,0,0,1,-3-h_c_f,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
t15=[0,0,0,0,0,0,0,0,0,1,0,0,0,1,-2-h_c_f,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
t16=[0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,-3,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0];
t17=[0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,-3-h_c_f,0,1,0,0,0,0,0,0,0,0,0,0,0,0];
t18=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,-3,1,1,0,0,0,0,0,0,0,0,0,0,0];
t19=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,-3-h_c_f,0,1,0,0,0,0,0,0,0,0,0,0];
t20=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,-3,1,1,0,0,0,0,0,0,0,0,0];
t21=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,-3-h_c_f,0,1,0,0,0,0,0,0,0,0];
t22=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,-3,1,0,0,0,1,0,0,0,0];
t23=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,-4-h_c_f,1,0,0,0,1,0,0,0];
t24=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,-3-h_c_f,1,0,0,0,1,0,0];
t25=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,-3-h_c_f,1,0,0,0,1,0];
t26=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,-2-h_c_f,0,0,0,0,1];
t27=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,-2-h_cc_f,1,0,0,0];
t28=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,-3-h_cc_f,1,0,0];
t29=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,-3-h_cc_f,1,0];
t30=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,-3-h_cc_f,1];
t31=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,-2-h_cc_f];

A=[t1;t2;t3;t4;t5;t6;t7;t8;t9;t10;t11;t12;t13;t14;t15;t16;t17;t18;t19;t20;t21;t22;t23;t24;t25;t26;t27;t28;t29;t30;t31];


B=[-h_a_f*Ta;
    -h_a_f*Ta;
    -h_a_f*Ta;
    -h_a_f*Ta;
    -h_a_f*Ta;
    0;
    0;
    0;
    0;
    0;
    0;
    -h_c_f*Tc;
    -h_c_f*Tc;
    -h_c_f*Tc;
    -h_c_f*Tc;
    0;
    -h_c_f*Tc;
    0;
    -h_c_f*Tc;
    0;
    -h_c_f*Tc
    0;
    -h_c_f*Tc;
    -h_c_f*Tc;
    -h_c_f*Tc;
    -h_c_f*Tc;
    -h_cc_f*Tcc;
    -h_cc_f*Tcc;
    -h_cc_f*Tcc;
    -h_cc_f*Tcc;
    -h_cc_f*Tcc];

X=linsolve(A,B)

a=rectangle('Position',[0 0 25 5]);
hold on
b=rectangle('Position',[0 5 5 20]);
c=rectangle('Position',[0 25 25 10]);
axis([-5 30 -5 40]);

Just need to plot it, either by using coordinates or by rectangle or any command

[KSSV]

If you have a matrix of temperature distribution, you can use pcolor to plot.