nonlinear coefficient in custom PDE.

I am trying to solve a variation of 2D steady state heat conduction PDE using PDE toolbox

Here , to simulate nonlinear behaviour.

Here is my code.

syms T(x,y)
syms k
pdeeq = diff(diff(T,x),x) + k*diff(diff(T,y),y)
symCoeffs = pdeCoefficients(pdeeq,T,'Symbolic',true)
symVars = [k]
symVals = [@(~,state) 0.01+1*state.u];
symCoeffs = subs(symCoeffs,symVars,symVals);

I get the following error, which does make sense

Error using sym/subs>normalize
Substitution expression X must be a symbolic, cell, or numeric array.

Error in sym/subs>mupadsubs (line 165)
[X2,Y2,symX,symY] = normalize(X,Y); %#ok

Error in sym/subs (line 153)
    G = mupadsubs(F,X,Y);

Error in sym/subs>@(value)subs(value,X,Y) (line 64)
        G = structfun(@(value) subs(value, X, Y), F, 'UniformOutput', false);

Error in sym/subs (line 64)
        G = structfun(@(value) subs(value, X, Y), F, 'UniformOutput', false);

In case of standard heat transfer I can write this.

k = @(~,state) 0.01+1*state.u
thermalProperties(model,'ThermalConductivity', k);

How to do the same with my custom PDE?