LJallo_f_q1 does not depend upon X(1), only upon X(2). If you happen to have the symbolic toolbox, you can find the expression for the results
If you assume that X(2) is non-negative, then the first element comes out (to one decimal place) with the form
-(2.7e+14*X2^(3/2))/((1.0e+15*X2 + 1.0)*(8.9e+3*X2^(1/2) + 1.0))
with non-negative X2, the numerator 2.7e+14*X2^(3/2) is positive and real-valued, and so is the denominator (1.0e+15*X2 + 1.0)*(8.9e+3*X2^(1/2) + 1.0) because each of the terms in the denominator is positive times X2 plus positive constant and that can never be less than the positive constant for non-negative X2.
So you have -positive/positive and that is -positive and you have to conclude that for any non-negative initial X2 that the function will continue to produce negative values that will decrease the boundary condition value, until eventually X(2) is negative.
Once you do not assume that X(2) is non-negative, then the first elemetn comes out (to one decimal place) with the form
-(1.2e+38*X2*conj(X2^(1/2)))/((5.4e+23*X2 + 5.2e+8)*(7.8e+18*conj(X2^(1/2)) + 8.8e+14))
Notice the X2^(1/2) . WIth negative X2 that is going to generate a complex result, and taking the complex conjugate of that is still going to be complex.
Putting these two together, any positive initial conditions are eventually going to lead to complex-valued results. The only question is whether your initial values are large enough that within your time span
If my calculations are right, then for that ic(1) value, your ic(2) would need to be no larger than about 6e-15 for you not to encounter complex values by time 1, but it is about 10^7 times larger than that.
Your plot does not show the plot starting from 0. Notice the scale of the plot on the y axes, that the marks are 1E-5 apart. Your plot starts from about 5E-8 which on the scale of the plot is not distinguishable from 0 unless you zoom in a fair way.
If you notice in the command window, after your plot you have
Warning: Imaginary parts of complex X and/or Y arguments ignored.
That is because your y is complex-valued.
By the way, please remember in your line of code
diffyq1=diffyq1';
that you are taking the conjugate transpose, not the regular transpose. That will not affect whether the results are complex or not, but it will affect the complex values returned. The plot of the imaginary component is considerably more interesting if you use the regular transpose .' instead of the conjugate transpose '
