theta34(:,i) = fsolve(@position, [1 1],options,th)*180/pi;
That th after the options struct will be interpreted according to a feature that has not been documented for quite a number of years. For a number of functions that do work on function handles, if you give extra parameters to the master function (fsolve(), fmincon(), ode45() and so on), then the extra parameters will be passed to the function handle.
So as well as the vector of the current theta value being passed to position(), the value of th will be passed as well.
The result is much like
theta34(:,i) = fsolve(@(theta)position)theta,th), [1 1],options)*180/pi;
except that the extra parameters are passed to all of the automatically invoked functions -- including any event functions or plotting functions or nonlinear constraints...
These days, if you want to pass an extra parameter to a function, use parameterization like I show here.
Meanwhile, your position function is not expecting th to be passed in.