save solution of integrant in 3d array

Hi, I got a problem and I'm relatively new to Matlab, so I'm not able to solve it myself. I'm working on the temperature distribution within a concrete sample caused by microwave power. This is a model i used to solve the heat transfer equation. with f_time(x,y,z) i'm hoping to get temperatures that i want to save in a 3d array, that i createdwith zeros over the length of my x,y,z coordinates. That does not work so far. Afterwards i want to slice plot different levels of the sample to compare it with my measured data. I'd be really thankfull for any help!This is the code i use:

syms X C;
a = 2;
c = 17; %seconds of microwave penetration
b = sqrt((1/4*c)+a);
m = (-1/(4*c))*(((3/4)-a*c)/(a*c+0.25));
q0 = 800; %in J/cm^2
T0 = 293.15; %Kelvin
xp= [0:40];
yp=[0:40]';
zp=[8:60];
[Xp, Yp, Zp] = ndgrid(xp, yp, zp);
f_time=zeros(length(xp),length(xp),length(zp));

for x=[0:3];
    for y=[0:3];
        for z=[8:11];
            X = [x;y;z];
            k = norm(X);
fun = @(C) ((a/(pi))^(3/2))*((1/(4*(pi)*C))^(3/2))*(b^(3/2))*exp((-m)*(k^2));
q =  integral(fun,0,c,'ArrayValued',true);
f_time(x+1,y+1,z) = T0 + 800*q;
save('f_time (x,y,z).mat');
        end
    end
end

[Mukul Rao]

Hi Julia, please provide more information on the equation you are trying to solve so that I can better relate to the problem at hand.