First, always test how your function operates. That will verify it produces what you expect.
Do you see it? theta returns a SCALAR RESULT.
Fsolve is a solver that tries to solve systems of equations. Here, you have 25 unknowns, as the elements of a 5x5 matrix. Can fsolve find a single set of parameters for that matrix, that will produce theta==0?
You have ONE equation that theta implicitly produces, but 25 unknowns.
Perhaps you want to use an optimizer, perhaps one that will minimize or maximize theta.
[xfinal,fval,exitflag] = fminunc(@theta,x0,opts)
So fminunc essentially found the function is unbounded. Apparently it could make that happen by setting x(5,5)=-1.
ffr = linspace(0.05,0.1,5);
omega = phi.*(1+psi)./(phi.*(1+psi)-1);
F = 1 + (x(i,j) - ffr(j))/(1+x(i,j))*omega(i)/(1-omega(i))+ psi(i)*((1+x(i,j))/(1+mean(x(i,j)))).^(omega(i)/(1-omega(i)));
Regardless, fsolve is arguably not the correct tool to solve your problem, since it cannot solve ONE equation with 25 unknowns. Since I don't really know what problem you wanted to solve, I am at a loss beyond this point.
I do have one serious issue with the function theta. Note that inside of the loop, you set F to a scalar value. It does so repeatedly. Perhaps you really intended to create a MATRIX F(i,j). So the ONLY value that ever gets returned is essentially F(5,5). And that explains why fminunc only bothered to iterate on the element x(5,5).
But this is just idle speculation on my part.