Two types of error when trying to solve time dependent pde

hello, i have written the following code, to solve a 3d pde.

the steady state solution is working properly, but fails to solve for the time dependent case.

whenever i type in the specifycoefficients line function handels, i get the following error: Solution does not correspond to time-dependent PDE.

and whenever i type values in this line (as shown in the attached code) i recieve this error: The value of 'colormapdata' is invalid. Operands to the and && operators must be convertible to logical scalar values. in this case my solution vector u does not converge to np*N size.

so i have 2 questions: why is there a difference between this two options? and what path should i take in order to fix this?

EDIT - i understand now that my solution vector u holds the solutions of every time lap. but i still don't understand why is there a difference when i type in the "specifycoefficients" line a value or a function handle.

thanks

d=.1; %depth
D=25*d/3; %diameter
ddr=d/D;  
a=D/2; %semi major axis

% generate alphashape
[az,el,r] = meshgrid(linspace(0,2*pi-.01,60),linspace(-pi,0,60),[0.99,1]);
[x,y,z] = sph2cart(az,el,r);
x=x*a;
y=y*a;
z=z*d+1;
shp = alphaShape(x(:),y(:),z(:),0.25);
% plot(shp);

%applying the geometry to the model
[elements,nodes] = boundaryFacets(shp);
nodes = nodes';
elements = elements';
N=1; %number of pde
model = createpde(N);
geometryFromMesh(model,nodes,elements);
% pdegplot(model,'FaceLabels','on','FaceAlpha',1);

% Generate the mesh
hmax=.025;
generateMesh(model,'hmax',hmax);
% pdeplot3D(model);

% Definition of PDE coefficients 
k = 5.e-3; % thermal conductivity, W/(m-degree C)
rho = 1500; % density, kg/m^3
cp = 500; % specific heat, W-s/(kg-degree C)

 cFunc = @(region,state) k*(region.x+region.y); 
  dFunc = @(region,state) rho*cp*(region.x+region.y);
  q = 0;
  fFunc = @(region,state) q*(region.x+region.y);
  A = 0;

% general properties and constants
sigSB = 5.670373e-8; % Stefan-Boltzmann constant, W/(m^2-K^4)
emiss = .95; % emissivity of the plate surfaces
lat=85; %defval('lat',85);
S0=1361; %solar constant W/m^2
albedo=0.1;
e=4; %sun elevation

%useful relations
w=tand(e);
p=d-D*w/2;
v=sqrt(-a^2*w^2+d^2);
s=w*a^2*p/v;
m=s^2+a^2*(d^2-p^2);

%fluxes
F=1/(1+(ddr^-2)/4); 
Fsolar=S0*(1-albedo)*cosd(lat);
Fingersoll=Fsolar*(F/(1-albedo*F))*(emiss+albedo*(1-F));

% Boundary Conditions
rad = @(region,state) emiss*sigSB*state.u.^3.*(region.x+region.y); 
solar = @(region,~) (region.x+region.y)* Fsolar;
shadow= @(region,~) (region.x+region.y)* Fsolar.*(((v*region.y+s).^2+d^2*region.x.^2)>m)+...
                    (region.x+region.y)* Fingersoll.*(((v*region.y+s).^2+d^2*region.x.^2)<m);

bbottom  = applyBoundaryCondition(model,'neumann','Face',2,'q',0,'g',0);
btop  = applyBoundaryCondition(model,'neumann','Face',1,'q',rad,'g',shadow);
bcircumference  = applyBoundaryCondition(model,'neumann','Edge',1,'q',rad,'g',solar);               

% Specify PDE Coefficients for Steady State Solution
% specifyCoefficients(model,'m',0,'d',0,'c',cFunc,'a',0,'f',0);

% try and solve the equation using solvepde

u0 = 100; %initial guess
setInitialConditions(model,u0);
% model.SolverOptions.ReportStatistics = 'on';
% results = solvepde(model);
% u = results.NodalSolution;

%solution with transient state 

% Specify PDE Coefficients for Transient Solution
specifyCoefficients(model,'m',0,'d',rho*cp,'c',k,'a',0,'f',0);

setInitialConditions(model,u0);
model.SolverOptions.ReportStatistics = 'on';
tfinal = 10000;
tlist = 0:100:tfinal;

result = solvepde(model,tlist);
u = result.NodalSolution;

% figure;
pdeplot3D(model,'ColorMapData',u);