Use solution of ODEs to fit experimental data: Invalid text character in solution

My code starts with solving a series of coupled differential equations, together with initial conditions, in order to find the constants for the solution. Of the solutions of all six equations, I focus on one of them p22(t), which is relatable to experimental data. To turn it into a correlation function, the thing that I measured, I have to divide the solution for p22(t) by the steady state solution p22(Inf).

syms p11(t) p22(t) p33xy(t) p33z(t) p12(t) p21(t) k21 k23xy k23z k31xy k31z K2 Gamma Omega 

%Defining the ordinary differential equations and load them into a matrix
ode1 = diff(p11) == k21*p22 + k31xy*p33xy + k31z*p33z - 1i*(Omega/2)*p12 + 1i*(Omega/2)*p21;
ode2 = diff(p22) == -(k21 + k23xy + k23z)*p22 + 1i*(Omega/2)*p12 - 1i*(Omega/2)*p21; 
ode3 = diff(p33xy) == k23xy*p22 - k31xy*p33xy; 
ode4 = diff(p33z) == k23z*p22 - k31z*p33z; 
ode5 = diff(p12) == -1i*(Omega/2)*p11 + 1i*(Omega/2)*p22 - Gamma*p12;
ode6 = diff(p21) == 1i*(Omega/2)*p11 - 1i*(Omega/2)*p22 - Gamma*p21; 
odes = [ode1;ode2;ode3;ode4;ode5;ode6];

% Making assumptions and setting initial conditions for the solution
assume(p11(t) + p22(t) + p33xy(t) + p33z(t) == 1); 
assume([k21 k23xy k23z k31xy k31z K2 Gamma Omega] > 0)

cond = [p11(0)==1, p22(0)==0, p33xy(0)==0, p33z(0)==0, p21(0)==0, p12(0)==0];

% Solving the system of ODEs
S = dsolve(odes,cond); 

% Solution of interest, to relate to the experimental autocorrelation
C = S.p22; 
C = C/(limit(C, t, Inf, 'left')); 

After execution I end up with a huge expression. However, substituting some numbers into the equation and plotting it shows me the desired result and behaviour while changing parameters.

The solution would probably have been even longer were it not for the countless expressions of the form:

root(#X^4 + (#X^3*(2*Gamma + 2*k21 + 2*k23z + 2*k31z + 2*k23xy + 2*k31xy))/2 + (#X^2*(2*Gamma*k21 + 2*Gamma*k23z + 2*Gamma*k31z + 2*Gamma*k23xy + 2*Gamma*k31xy + 2*k21*k31z + 2*k21*k31xy + 2*k23z*k31z + 2*k23z*k31xy + 2*k31z*k23xy + 2*k31z*k31xy + 2*k23xy*k31xy + 2*Omega^2))/2 + (#X*(Omega^2*k23z + 2*Omega^2*k31z + Omega^2*k23xy + 2*Omega^2*k31xy + 2*Gamma*k21*k31z + 2*Gamma*k21*k31xy + 2*Gamma*k23z*k31z + 2*Gamma*k23z*k31xy + 2*Gamma*k31z*k23xy + 2*Gamma*k31z*k31xy + 2*Gamma*k23xy*k31xy + 2*k21*k31z*k31xy + 2*k23z*k31z*k31xy + 2*k31z*k23xy*k31xy))/2 + Gamma*k31z*k23xy*k31xy + Gamma*k23z*k31z*k31xy + Gamma*k21*k31z*k31xy + (Omega^2*k31z*k23xy)/2 + (Omega^2*k23z*k31xy)/2 + Omega^2*k31z*k31xy, #X, 3)

When I want to fit formula C to my experimental data, I tried to use the Curve fitting tool. It gave me the Error: Invalid text character, maybe related to the hashtag. Also doing it in the code by trying to convert C to a matlabFunction I get an error saying: Code generation failed due to unexpected object of type 'RootOf'. Why is there a #X in the expression and how can I fit my formula to some data?