You cannot provide an inequality as a bounds for int(). You can do indefinite integration by not providing a bound at all, and you can do definite integration by providing a lower bound and an upper bound, but you cannot provide an inequality and have MATLAB solve it for the domain of validity and integrate over that domain. (For one thing, defining a domain of validity does not define an integration direction, does not define which is the upper bound and which is the lower bound.)
What you need to do is to solve() the inequality first, and extract the lower bound and upper bound of range the solution is indicated as valid for, and use those bounds for your integration. See the 'ReturnConditions' option of solve() in regards to that. You might find that in practice you need to solve() as an equality rather than as an inequality, and then figure out whether the solution is the upper or lower bound and use it appropriately in integration.
Your L1 comes out as -1/(-1+erf((2/3)*y*sqrt(2))) and your L2 comes out as 1/(erf((2/3)*y*sqrt(2))+1) . If you solve(L1=T,y) then the value returned is the lower bound for L1>T. If you solve(L2==0,T) then the value returned is the lower bound for L2<T.
I notice that your forms are exactly the same for error1 and error2, with the only difference being which side of the inequality you are choosing. I recommend that you merge both of them into the same routine that returns two values, instead of having to go through all the hard work twice.
int(I2,y) has a nice closed form representation, but int(I1,y) does not, so your alpha and beta are going to be represented by expressions that include unevaluated integrals even though you may have given definite T values. That basically indicates that a numeric integral must be done at that point.
What I suggest you do is calculate error1 and error2 on a symbolic T (and merge the two into a single routine), returning symbolic formula. Then use matlabFunction() to convert the symbolic formula into an anonymous function suitable for being invoked with a definite numeric T and returning a numeric (not symbolic) result.
matlabFunction() normally produces vectorized code, so normally you would be able to invoke it on your t = 0:0.01:30 and get back a vector of results that you could plot. I am, however, not certain that matlabFunction will be able to vectorize when numeric integration is involved. I do not have the Symbolic Toolkit so I cannot test that myself. If not then the adjustments are not difficult.
Caution: the solution to the inequalities are only valid as upper or lower bounds for the case of T >= 1/2. Your L1 and L2 both have lower bounds of 1/2 for real-valued y. You need to decide whether you want to be integrating over complex y or not -- nothing in your problem statement disallows doing so, but when your boundaries are complex valued, you cannot uniquely define upper and lower bounds for the inequalities and so cannot uniquely say what you are integrating over.
Best Answer