Your fun3 line includes
symsum(S0i(ik)*fun0(1i*PSIi{ik}-(Gamma1-1i*delta)*SPSI), ik, 2, n)
Inside the symsum, you try to use the symbol ik to index PSIi .
In MATLAB, it is never possible to use a symbol to index something -- not possible for vectors, not possible for cell arrays.
I have attached a re-worked version. It is not possible to do the integration using quadqk because at the point of the integration you still have an unknown parameter, THETA1, that you will be one of the two parameters you are minimizing over.
I would note, though, that your exp(-(delta*x1)) and your max(q,0) are separable, so you can reduce the work required for minimization. exp(-(delta*x1)) is strictly >= 0 for real-valued x1. max(q,0) is >= 0 for all q. The product must then be >= 0. You could achieve the minimum of the product if either value were 0.
q looks complicated to compute and solve for a root, but the root of exp(-(delta*x1)) is easy to compute: for positive delta, it occurs when x1 is infinity, giving exp(-infinity) which is 0. Your fminbnd includes the entire range -inf to +inf so +inf is a valid value for x1.
Therefore, by examination we can solve the entire problem as x = +infinity, without ever having computed q .
The only exception would be if we could prove that q is everywhere +infinity, over all THETA1, as in that case you would be working with zero times infinity which is undefined. Therefore if we could prove that q is +infinity for all THETA1, we would have to prohibit x1 = infinity, and so exp(-(delta*x1)) would be positive; in that case, since the positive value is being multiplied by infinity (by hypothesis), it does not matter what value of x1 we use (other that +infinity) as all of them would minimize the product to the same value, +infinity.
You should probably proceed from here to prove that at least one value of q is finite for some THETA1, to make the overall solution x1 = +infinity valid.
