matlabFunction returns a bad function

I am using Symbolic toolbox in order to use ztrans & iztrans.

Here's the code for computation, where P is some diagonalizable matrix:

function Tij=getTdistribution(P, i, j)
[V, D] = eig(P);
syms n
syms Dn
Vinv = inv(V);
Dn = sym(D)^n;
P = V*Dn*inv(V);
pij = P(i, j);
pjj = P(j, j);
Tij = iztrans(ztrans(pij)/ztrans(pjj));
end

The problem is that when I call

>> Tij = getTdistribution(P, i, j)
>> f = matlabFunction(Tij)

and try to evaluate f (e.g. f(0), f(1)), I get the following error:

>> f(1)
Unrecognized function or variable 'r6'.

This is weird because I am expecting f to only depend on n.

Please explain why this is happening and how to fix this :)

[Stephan]

With which P, i and j do you call the function? For me this example worked in R2020a:

P = [4 -3 -3; 3 -2 -3; -1 1 2]
ii = 1
jj = 2

Tij = getTdistribution(P, ii, jj)
f = matlabFunction(Tij)
f(1)


function Tij=getTdistribution(P, ii, jj)
[V, D] = eig(P);
syms n
syms Dn
Vinv = inv(V);
Dn = sym(D)^n;
P = V*Dn*Vinv;
pij = P(ii, jj);
pjj = P(jj, jj);
Tij = iztrans(ztrans(pij)/ztrans(pjj));
end

result is:

P =

     4    -3    -3
     3    -2    -3
    -1     1     2


ii =

     1


jj =

     2

 
Tij =
 
((287115466509349393797924020813824*2^(1/2)*7^(1/2) - 1074287706121754825198050084043197)*kroneckerDelta(n, 0))/(430673199764024090696886031220736*14^(1/2) - 1327477824731508952975575249674352) + ((((17944716656834337112370251300864*14^(1/2) - 37564467627284328086361790527834)/(8972358328417168556185125650432*14^(1/2) - 27655788015239769853657817701549))^n/(17944716656834337112370251300864*14^(1/2) - 37564467627284328086361790527834) - kroneckerDelta(n, 0)/(17944716656834337112370251300864*14^(1/2) - 37564467627284328086361790527834))*(318466831773461748772018933149303065732501697143136955815428096*2^(1/2)*7^(1/2) - 477700247660192623158028399723954598598752545714705433723142144*14^(1/2) + 280837030099072488269234933106573354835797106198448123800453120))/(26917074985251505668555376951296*14^(1/2) - 82967364045719309560973453104647)
 

f =

  function_handle with value:

    @(n)((sqrt(2.0).*sqrt(7.0).*2.871154665093494e+32-1.074287706121755e+33).*(n==0.0))./(sqrt(1.4e+1).*4.306731997640241e+32-1.327477824731509e+33)+((((sqrt(1.4e+1).*1.794471665683434e+31-3.756446762728433e+31)./(sqrt(1.4e+1).*8.972358328417169e+30-2.765578801523977e+31)).^n./(sqrt(1.4e+1).*1.794471665683434e+31-3.756446762728433e+31)-(n==0.0)./(sqrt(1.4e+1).*1.794471665683434e+31-3.756446762728433e+31)).*(sqrt(2.0).*sqrt(7.0).*3.184668317734617e+62-sqrt(1.4e+1).*4.777002476601926e+62+2.808370300990725e+62))./(sqrt(1.4e+1).*2.691707498525151e+31-8.296736404571931e+31)


ans =

   -3.0000