Adding heat source term in LBM

Okay, I looked up almost every where and i could not find any answer to this question.
(1) you have 1 dimensional diffusion equation that you want to solve with LBM.
If say, you choose to do it with D1q2 or d1q3 for the sake of simplicity.
My west and east boundaries are at 100 and 200 Celsius.
A.A mohammed book says that you have to add this term to the source code Q_generated/(rho*Cp);
k is given as 0.5 w/mk.
Q_gen=1000 KW.
Thickness =2 cm.
for the love of god,can someone please solve this for me. i cannot get the solution for this simple problem.
I’ll give here the analytical solution with 5 inner nodes
T profile=[100 146 214 250 254 226 200]
Many papers suggest that adding heat source just adding a number multiplied by weighing factor. But its simply not working.
My matalb code:

clear all;
close all;
clc;

% Intial values
tw=1;
te=0;
w=[0.5 0.5];
n=100;
% g(:,:,1)=w(1)teones(1,6);
% g(:,:,2)=w(2)teones(1,6);
alpha=1;
omegat=1/(alpha+0.5);

T(1)=100;
T(n)=200;
T(2:n-1)=0;
tw=T(1);
te=T(n);

 %equillibrium

%
for i=1:2
g_eq(:,:,i)=w(i)*T;
end
g=g_eq;
tsteps=4000;

qg=20000;
% cp=100;
source=0;
for ii=1:tsteps
ii

 %collision
 rho=sum(g,3);
 
 
 for i=1:2
     % here is where i add the source term. how do i calculate this
     % number? 
     g(:,:,i)=g(:,:,i)-omegat*(g(:,:,i)-g_eq(:,:,i))+w(i)*0.9;
 end
 
  %stream
 g(:,:,1)=[g(:,1,1) g(:,1:n-1)];
 g(:,:,2)=[g(:,2:n,2) g(:,n,2)];
 

 %BC
 g(:,1,1)=w(1)*tw+w(2)*tw-g(:,1,2);
 g(:,n,2)=w(1)*te+w(2)*te-g(:,n,1);
  T(:,:,ii)=sum(g,3);

% source=qg./rho;
if ii>1
delta(:,:,ii)=T(:,:,ii)-T(:,:,ii-1);
if (delta(:,:,ii)) < 10^-6
break;
end
end

% T=sum(g,3);

 for i=1:2
 g_eq(:,:,i)=w(i)*T(:,:,ii);
 end
plot( T(:,:,ii),'*')

end

% axis([1 n tw te])