You can solve diffusion equation in PDE Toolbox. Check this example:
% Create a model.
model = createpde(1);
% Create geometry and assign it to the model
R1 = [3;4;-1;1;1;-1;-1;-1;1;1];
C1 = [1;0;0;0.4];
C1 = [C1;zeros(length(R1) - length(C1),1)];
gd = [R1,C1];
sf = 'R1+C1';
ns = char('R1','C1')';
g = decsg(gd,sf,ns);
geometryFromEdges(model,g);
figure
pdegplot(model,'EdgeLabels','on','FaceLabels','on')
% Specify coefficients, diffusion coefficient as c and source as f
% Assign zero source and a diffusion coefficient, c = 1
specifyCoefficients(model,'Face',1,'m',0,'d',1,'c',1,'a',0,'f',0);
% Assign non-zero source and a different diffusion coefficient c = 2 on inner
% circle
specifyCoefficients(model,'Face',2,'m',0,'d',1,'c',2,'a',0,'f',1);
% Set initial condition
setInitialConditions(model,0);
setInitialConditions(model,1,'Face',2);
% Generate mesh and solve.
generateMesh(model)
tlist = linspace(0,0.1,20);
result = solvepde(model,tlist)
% Plot results of last time-step.
figure
pdeplot(model,'XYData',result.NodalSolution(:,end))
