Joule Heating in a single cell

Hello. I am having a problem with my code. I have a MATLAB code for joule heating in a single cell isoelectric focusing and I want to measure the temperature changes with electrodes close

to each other and far from each other. Currently in this model I cant really see termperature difference but in practical model temperature rises to about 5 degree. I am stuck at this point in coding.

The rate of heat generated (Q) due to a resistor of length L and cross sectional area A, with an elementary resistivity of γ is given by,

(1)

Where E is the applied electric field, and V is the applied potential. Heat is dissipated from the device throughout electrophoresis. The rate of heat dissipation via conduction (P) can be generally expressed as,

(2)

Here T is the higher temperature in the device that drives heat conduction, Tw is the ambient temperature, κ is the thermal conductivity of the material, L2 is the distance of heat dissipation (here it is proportional of the thickness of the device), and S is the heat dissipating area, (here it is the top / bottom surface area).

I equated the above two equations and plotted y = mx+c

I have attached my code. Can someone help me in this matter. Thanks in advance.

% x = 0 : 0.001: 1.2 ; %X is the length 40 Microns of Gel thickness plus half a milimiter for Lid
% Polyacrymide =  1.2 : 0.001: 1.24;
% Agarose = 1.24 : 0.001: 1.34;
% Glass = 1.34 : 0.001: 2.54 ;
% 
% E = 50; % 50V/CM Applied electric field
% V = 150; % 150V Applied potential over 3CM
% 
% Length = 75; %75mm
% Width = 25; %25mm
% Height = 0.140; %0.140 micrometer
% 
% S = (2*Length*Width) + (2*Length*Height) + (2*Height*Width); %Heat dissipating area
% K = 0.56; % W m −1 K −1 Thermal conductivity of material https://www.researchgate.net/publication/10593278_Measurement_of_the_thermal_conductivity_of_polyacrylamide_tissue-equivalent_material
% A = 75*25; % Cross-sectional area 
% G = 6; % Elemental resistivity
% m1 = E*V*A; 
% m2 = K*S*G;
% m = m1/m2; % Slope
% 
% Tw = 293; %Tw is the ambient temperature and is the intercept; Kelvin
% T = m*x + Tw; %T is the Y axis in this instance
% 
% figure(1)
% plot(x, T)
% % grid
% % axis([0  5    0  12])

%should come out to 2K 
L = 140*10^-6; %m
E = 50*10^2; %V/m
V = 150; %V
A = 140*10^-6*25*10^-3; %m^2
k = 1.46; %W/(m*K)
S = 2*(75/2)*25*10^-6; %m^2 
%electrical resistivity equation
gamma = (150)*(A)./((10*10^-3)*(3*10^-2)) %ohm m
% gamma = 1/.09; 

L*E*V*A./(k*S*gamma)

%should come out to 90K 
L = ((5*10^-3)+(40*10^-6)); %m
E = 50*10^2; %V/m
V = 150; %V
A = L*25*10^-3; %m^2
k = 1.46; %W/(m*K)
S = (75/2)*25*10^-6; %m^2 
%electrical resistivity equation
gamma = (150)*(A)./((250*10^-3)*(3*10^-2)) %ohm m 
% gamma = 1/.09; %Julea paper

temperature_change = L*E*V*A./(k*S*gamma)

%scIEF

V = 690;
I = 0.22*10^-3;
A = 25*10^-3*0.5*10^-3;
L = 3*10^-2;
gamma = V*A./(I*L)

L = .5*10^-3; %m
E = 690/(3*10^-2); %V/m
V = 690; %V
A = L*25*10^-3; %m^2
k = 1.46; %W/(m*K)
S = 10*25*10^-6; %m^2 

L*E*V*A./(k*S*gamma)

room_temp = 293
gel_temp = 293 + L*E*V*A./(k*S*gamma)

Output:-

>> Simulation_gl

gamma =

1.7500

ans =

0.0767

gamma =

2.5200

temperature_change =

138.0822

gamma =

1.3068e+03

ans =

0.2079

room_temp =

293

gel_temp =

293.2079

As you can see, room_temp and gel_temp are almost same but in practical there is about 5 degrees of difference.