Sol_p=solve(ycl_p-ygl_p,xint2); is going to have two results. You then assign it to xint_cg, so xint_cg is going to have length 2. But your "for Il" loop tries to index it up to numel(k) which is 49001. That would fail as soon as I1 became 3.
In the line
SOl(Il)=solve(e1(Il)-e2(Il),xint_cg(Il));
xint_cg(I1) is going to be an expression, and e1(Il)-e2(Il) is going to be an expression, so the solve() will be trying to solve two expressions simultaneously. Which two variables are you expecting to solve over? Even if the expressions do happen to have exactly two free variables there, the code would be so much more understandable if you named the specific variables there. Do not make life harder than it needs to be for people who are reading your code.
You are trying to find the Ilmax that makes xs 0, but nothing in your code sets I1max or uses it (other than that one line), so we do not know what you are looking for. Do you perhaps want to assign Il there? Are you looking for the Il that makes xs 0? Is there reason to expect that the solution will be an integer?
Why are you defining
k=1000:1:50000;
and iterating Il from 1 to numel(k), but your only other reference to k is to assign Ilmax to it? Is your code possibly needing to refer to k(Il) in those places it is referring directly to Il (other than as a subscript) ?
When you do assign the (presumably scalar) Ilmax to k, are you expecting that the "for Il=1:numel(k) will notice the change to the size of k and so exit the loop, numel(k) having become 1? Because that is not how "for" loops work: for loops record copies of the loop parameters and use those copies no matter how any variables might change in the loop.
Note: it is not possible to solve() with respect to an expresion, only with respect to a variable.
