vpasolve for nonlinear system of equations instead of fsolve for complex valued results

I have nonlinear 2 quite complex equations for 2 unknowns. At first, I used fsolve. But, because the equations are complex valued, I was taught vpasolve should be used instead of fsolve.

The following is the code for this problem.

X = sym('x', [1,2]);
T = 1/100;
F = torder1(X,T);
sols = vpasolve(F, X);
sols.x1
sols.x2

Below is torder1.m function called from above script.

function F=torder1(x,T)
   x_1=[0:0.01:1];
   b=0.6; 

syms x_1 h theta 

f_11 = 1-( (h+(x_1-b)*theta)^2/(h+(x_1-b)*theta-1*x_1*(1-x_1))^2 );
f_21 = (x_1-b)/2*( 1-( (h+(1-b)*theta)^2/(h+(x_1-b)*theta-x_1*(1-x_1))^2 ));

fint_1 = int(f_11, x_1);
fint_2 = int(f_21, x_1);

x_1=1;
upper_1=subs(fint_1);
upper_2=subs(fint_2);
clear x_1;
x_1=0;
lower_1=subs(fint_1);
lower_2=subs(fint_2);
clear x_1;
integral_result_1old=upper_1-lower_1;
integral_result_2old=upper_2-lower_2;

h0 = kappa *b*(1-b);
theta0 = kappa*(1-2*b);

integral_result_1 = subs(integral_result_1old, {h, theta}, {x(1), x(2)});
integral_result_2 = subs(integral_result_2old, {h, theta}, {x(1), x(2)});
% 
% expression1 = double(x(1) - integral_result_1*T^2 -h0);
% expression2 = double(x(2) - integral_result_2*T^2 - theta0);

F = [x(1) - integral_result_1*T^2 -h0;
    x(2) - integral_result_2*T^2 - theta0];

% F = [double(x(1) - integral_result_1*T^2 -h0);
%     double(x(2) - integral_result_2*T^2 - theta0)];

The errors are two types:

If I use F with double inputs inside torder1.m, when I run the script I get

DOUBLE cannot convert the input expression into a double array. If the input expression contains a symbolic variable, use VPA.

If I use F without double arguments inside torder1.m, I get :

CAT arguments dimensions are not consistent. error.

You might would like to see earlier question related to this problem

http://uk.mathworks.com/matlabcentral/answers/300450-matlab-limitation-in-fsolve-using-function-input#answer_232497

Looking forward to hearing any help.. Thanks.