There are smart ways of approaching this problem, and there are silly ways of approaching this problem.
One option of getting around Matlab imaginary joint values is to implement the code in python with regular floating point values and then import the package into matlab. That way, Python will throw errors when complex joint values are returned (rather than propagate complex values through the solver).
A more correct solution to your problem would be to implement datatype checking all throughout the solver. Each step of the Ik solver should check if you have imaginary values. Then, once you have found the imaginary value, you can figure out where the imaginary values are being generated. The most common problem in my experience comes from the use of inverse trig functions being given values outside of their domain. When I get imaginary values, I am passing a value of magnitude great than 1 into an inverse trigonometric functions like acos(). Matlab handles this problem by returning complex values, and that does not help in debugging. That is why I recommend moving to python (or better yet, C++). You might also be running into division by 0, but the most likely issue is that you are calling trig and inverse trig functions with improper inputs.