Solving multiple variable functions: how to create a symbolic array?

I am using the following functions to model a reactor of fixed volume, with inlet temp Ti, operating temp To and proportion of reactant in feed y_A. The reactor has an operating point dictated by the energy balance dH_gen=dH_rem. Over a range of To and y_A I want to solve dH_gen-dH_rem==o to find Ti at the operating point using the follow matlab code.

%defining constants
V=8.5;
P=3e5; %Pa
R=8.3145; %kJ/mol.K
Cp=1.09 %kJ/kgK
Mr_air=29e-3; %kg/mol
dH_1=-1850; %kJ/mol
dH_2=-1423;
dH_3=-3273;
F_ao=0.1; %molar flow rate of benzene

Ti=600:1:850; 
y_A=0.001:0.0001:0.01;
hold on;
for i=1:numel(Ti)

      for j=1:numel(y_A)
       To = sym('To', [i j]);
      k1(i) = 1e7.*exp(-12700./To);
      k2(i) = 5e4.*exp(-10800./To);
      k3(i) = 7e7.*exp(-15000./To); 
          t(j)=(V.*y_A(j).*P)./(F_ao.*R.*To);
          X(j)=(t(j).*(k1(i)+k3(i)))./(1+t(j).*(k1(i)+k3(i)));
          Y_B(j) = (k1(i).*t(j))./(((k2(i).*t(j))+1).*(1+(t(j).*(k1(i)+k3(i)))))
          dH_gen(i,j)=-((V.*y_A(j)/F_ao.*P)./(R.*To)).*((((k1(i).*dH_1)+(k3(i).*dH_3)).*(1-X(j)))+(k2(i).*Y_B(j).*dH_2)); %enthalphy generated by exothermic reaction
          dH_rem(i,j)=(Mr_air./y_A(j)).*Cp.*(To-Ti(i));%enthalphy change of process fluid
          eqn(i,j) = dH_gen(i,j)-dH_rem(i,j);
          solx = fzero(eqn(i,j),To)
      end
  end

From this, there is an error occurring when converting from sym to double: DOUBLE cannot convert the input expression into a double array. What does this mean and how do I fix this?