sol_tau =
I think this is what you want. But although there are four boundary conditions for four free parameters, they do not suffice to fix a unique solution. Thus there must be an additional condition missing.
syms To(x) Ti(x) tau(x)
syms beta c tau_avg
% Define ODE system and boundary conditions
eqns = [diff(To,x)+tau==0,diff(Ti,x)-2*tau==0,diff(tau,x,2)-beta^2*tau==0];
conds = [To(-c) == 2*c*tau_avg,Ti(-c) == 0, To(c) == 0, Ti(c) == 4*c*tau_avg];
% Compute solution of ODE system with boundary conditions
sol = dsolve(eqns,conds);
sol_tau = simplify(sol.tau)
% Write solution for tau in basis functions sinh(beta*x) and cosh(beta*x)
syms A B
eqn = sol_tau == A*sinh(beta*x)+B*cosh(beta*x);
coeffs = solve([subs(eqn,x,c),subs(eqn,x,-c)],[A,B]);
coeffs.A = simplify(coeffs.A);
coeffs.B = simplify(coeffs.B);
sol_tau = coeffs.A*sinh(beta*x) + coeffs.B*cosh(beta*x)
v4 = symvar(sol_tau)
% Choose numerical values for beta, tau_avg, c and the free paramter C1 and substitute
% in solution for tau
beta_num = sqrt(0.048);
tau_avg_num = 0.25;
c_num = 1;
C1_num = 0;
sol_tau_num = subs(sol_tau,[beta tau_avg c v4(1)],[beta_num tau_avg_num c_num C1_num])
% Plot solution for tau
fplot(sol_tau_num,[-c_num c_num]) % plot tau_c




