Two dimensional transient heat conduction problem with heat source

Question

There is a question about 2D heat conduction with heat source in the center. and its image is like this. Given: The 2cm*2cm heat source could be considered as one node. The temperature of all the nodes at t=0s is know, and the generation heat flux q_generation is also knowed. None of the numbers in the question are important. Find: the temperature distribution of the 2D plane changing with time.

I have known that the numerical method of 2D(node by node) without heat source already, but if there is a heat source which is generating heat, I don't know how to solve this problem. where should I start to calculate? I make a MATLAB code, and get the apparently wrong answer.

format longG;
clc;clear;
%% Define Parameters and Initial Conditions
Lx = 0.3; % Length of the region
Ly = 0.2; % Width of the region
Nx = 31; % Number of points in each row/column
Ny = 21; % Number of points in each column/row
Tmax = 1; % Total time
dt = 10^(-2); % Time step
num_time_steps = Tmax/dt; % Total number of time steps
k = 14.4; % Thermal conductivity
h = 100; % Heat transfer coefficient
alpha = 0.387*10^(-5); % Thermal diffusivity
Cp = 461; % Specific heat capacity
den = 7817; % Density
q_gen = 50*10^(-6)*0.01*1;  % Heat source (50 watts per cubic centimeter)
num_i = 0;
num_j = 2;

dx = Lx / (Nx-1);
dy = Ly / (Ny-1);

%% Create an array to store temperature distribution at all time points
T = zeros(Ny, Nx, num_time_steps);
% Initial temperature distribution at t=0 (k=1)
T(:, :, 1) = 25*ones(Ny, Nx);

%% Time Stepping
% Loop to calculate temperature distribution at each time step
for t = 2:num_time_steps
    for i = (Ny+1)/2:Ny  % Horizontal from 11-21
        for j = (Nx+1)/2:Nx % Vertical from 16-31
            % Use the previous time step array to compute the next time step array
            if i ~= 11
                num_j = 2;
            end
            if i == (Ny+1)/2 && j == (Nx+1)/2
                % Calculate temperature at the corner
                T(i, j, t) = ((T(i-1, j, t-1)-2*T(i,j,t-1)+T(i+1,j,t-1))/dx^2 + (T(i, j+1, t-1)-2*T(i,j,t-1)+T(i,j-1,t-1))/dy^2 + q_gen/k)*alpha*dt+T(i,j,t-1);
            elseif i == Ny && j == Nx % Corner
                T(i, j, t) = ((0.5*k/dy^2+0.5*k/dx^2+h/(dy*2)+h/(dx*2))*-T(i, j, t-1)+...
                    (0.5*k/dy^2*T(i,j-1,t-1)+0.5*k/dx^2*T(i-1,j,t-1)+h/(dy*2)*25+...
                    h/(dx*2)*25))*dt+T(i,j,t-1);
                T(i-num_i,j,t)=T(i,j,t);
                T(i,j-num_j,t)=T(i,j,t);
                T(i-num_i,j-num_j,t)=T(i,j,t);
                message = sprintf('Corner has been processed, %d, %d', i, j);
                disp(message);
            elseif i == Ny && j ~= Nx % Rightmost column
                T(i, j, t) = ((alpha/dx^2+0.5*alpha/dy^2+0.5*alpha/dy^2+h/(dx*den*Cp))*-T(i, j, t-1)+...
                    (alpha/dx^2*T(i-1,j,t-1)+0.5*alpha/dy^2*T(i,j+1,t-1)+0.5*alpha/dy^2*T(i,j-1,t-1)+...
                    h/(dx*den*Cp)*25))*dt+T(i,j,t-1); 
                T(i,j-num_j,t)=T(i,j,t);
                message = sprintf('Right column has been processed, %d, %d', i, j);
                disp(message);
            elseif j == Nx && i ~= Ny % Bottommost row
                T(i, j, t) = ((alpha/dx^2+alpha/dy^2+h/(dy*den*Cp))*-T(i, j, t-1)+...
                    (0.5*alpha/dx^2*T(i-1,j,t-1)+0.5*alpha/dx^2*T(i+1,j,t-1)+alpha/dy^2*T(i,j-1,t-1)+...
                    h/(dy*den*Cp)*25))*dt+T(i,j,t-1);
                T(i-num_i, j, t) = T(i,j,t);
                message = sprintf('Bottom row has been processed, %d, %d', i, j);
                disp(message);
            else
                T(i, j, t) = ((T(i-1, j, t-1)-2*T(i,j,t-1)+T(i+1,j,t-1))/dx^2 + (T(i, j+1, t-1)-2*T(i,j,t-1)+T(i,j-1,t-1))/dy^2)*alpha*dt+T(i,j,t-1); % Update temperature with arbitrary values
                if i == 11
                    T(i,j-num_j,t)=T(i, j, t);
                    num_j = num_j + 2;
                else
                    T(i-num_i,j,t)=T(i, j, t);
                    T(i-num_i,j-num_j,t)=T(i,j,t);
                    T(i,j-num_j,t)=T(i,j,t);
                    num_j = num_j + 2;
                end
                message = sprintf('Inner region has been processed, %d, %d', i, j);
                disp(message);
            end
        end % Calculation for each row
        num_i = num_i + 2;
    end % Calculation for one time step

% Plot Temperature Distribution
    % Create a figure window
    figure(1);
    imagesc(T(:,:,t));
    colorbar;
    title(['Time = ' num2str(t*dt)]);
    xlabel('X');
    ylabel('Y');
    drawnow;

    num_i = 0;
    num_j = 2;
end