Here's a brute force approach that at least seems to give a reasonable result for V0 = 0. This code runs very slowly because erfi() with a double input still is evaluated via the Symbolic version of erfi().
V0 = 0;
% Calculated
%
%% Solving the Maier-Saupe (MS) equation for the order parameter S
% Setting up the equations
syms x OP T_star
numerator = int((3*x^2/2-1/2)*exp((5*OP-V0)*(3*x^2/2-1/2)./T_star), x, 0 ,1);
denominator = int(exp((5*OP-V0)*(3*x^2/2-1/2)./T_star), x, 0, 1);
MS = (numerator./denominator) - OP;
func = matlabFunction(MS);
% Solving OP for each value of T*
OPval = -1:1e-2:1;
MS as defined above has a singularity at OP = 0, so remove OP = 0 from the values to be evaluated. Note that this will cause a root to be found at OP = 0, don't know if that's correct or not.
OPval(OPval == 0) = [];
T_star_min = 0.1;
T_star_max = 2.0;
interval = 0.01;
T_star_val = T_star_min:interval:T_star_max;
S = nan(numel(T_star_val),3);
for ii = 1:numel(T_star_val)
y1 = func(OPval,T_star_val(ii));
y = [0; y1(:)];
ys= [y1(:); 0];
k = find(sign(y) .* sign(ys) == -1);
OPzeros = -y1(k-1).*(OPval(k)-OPval(k-1))./(y1(k)-y1(k-1))+OPval(k-1);
S(ii,1:numel(OPzeros)) = sort(OPzeros);
end
S
plot(T_star_val,S,'bo')
