Thermal distribution problem in a iron bar in Matlab

I have to solve a 2D Steady State Heat Conduction Problem with heat generation (in Cartesian Coordinates) using Finite Difference Method.

I'm trying to draw isothermal lines and heat flow lines with Matlab.

Here is the link for one similar problem in 2D steady state heat transfer, but without heat generation: https://octave-online.net/bucket~cSmKVS6hKCGAPBW26Bj4S

In my case I have a square cross section bar 0.02m in the interior of which we have heat production, q=50 MW/m3, from a square section B, each side of which is H=L/2=0.01m.

We know that dx=dy=0.001m, distance between the nodes. For dx=dy my equation of Finite Difference Method is :

T(i-1,j) + T(i+1,j) - 4T(i,j) + T(i,j-1) + T(i,j+1) = -(q/k)*dx^2;

Outer surface temperature: T=300K=constant

Thermal conductivity: k=40W/mK

Density: ρ=7500kg/m3

Coefficient of thermal conductivity: cp=125J/kgK

Here is an image of the problem:

Here is my attempt, this code is based on the link above with modifications in SteadyState2D.m :

clear;clc;hold off;close all
q=50; 
k = 40; 
Lx = 0.02; 
Ly = 0.02; 

T_top    = 300 ; 
T_bottom = 300 ; 
T_left   = 300 ; 
T_right  = 300 ;  
dh=0.001; 
 

x=0:dh:Lx; x(end)=Lx;
y=0:dh:Ly; y(end)=Ly;
nx=length(x);  
ny=length(y);  
n=nx*ny;  

figure(1);plotMesh(x,y)
 
top=ny:ny:n;
bottom=1:ny:nx*ny;
left=1:ny;
right=(nx-1)*ny+1:n;

inner=1:n;
inner([left right top bottom])=[];
 
A=sparse(n,n); 
b=zeros(n,1);  

for i=inner
    A(i,i-ny)=1;
    A(i,i-1)=1;
    A(i,i)=-4;
    A(i,i+1)=1;
    A(i,i+ny)=1;
    
    b(i)=-q/k*(dh)^2;
end

for i=top
    A(i,:)=0; A(i,i)=1;
    b(i)=T_top;
end

for i=bottom
    A(i,:)=0; A(i,i)=1;
    b(i)=T_bottom;
end

for i=left
    A(i,:)=0; A(i,i)=1;
    b(i)=T_left;
end
for i=right
    A(i,:)=0; A(i,i)=1;
    b(i)=T_right;
end