You CANNOT fit that model. Why not? a and b are not independent parameters. They look as if they are.
Ignoring the part where x is small, your model looks like:
a*b*sqrt(x^2 - 1)
First, I'll point out that it does not even fit your data. When x is near 1, your model predicts ZERO. Yet your data passes through 1 or -1 when x is near zero. So your model does not even have the correct behavior for x near 1. Look at the curve you want to use. a and b are irrelevant, as they just scale the thing up and down.
As you should see, your data is flipped around. While both curves seems to be roughly asymptotically linear for large x, your data has a roughly zero slope near x == 1, whereas the model you want to use has an infinite slope near x==1.
But disregarding the fact that your model is complete crap for this data, could you even solve for both a and b?
Suppose you did find an optimal value for a and b. Now, consider what would happen if you tried instead a*2 and b/2 in that model? You would then notice that:
a*b*sqrt(x^2-1) == (a*2)*(b/2)*sqrt(x^2-1)
so if you think you have found an "optimal" set of values for a and b, note that they are not optimal. There are infinitely many possible sets of numbers you could have chosen, all equally optimal.
And that means you can arbitrarily set b==1. Now, in theory, you could hope to estimate a.
Except of course, that this model will never look like that data. lsqcurvefit will not help to fit a square peg into a round hole.
(Sigh. Do I need to tell you how to fix this? Probably.)
You seem to have the x and y relationships swapped around. That is, suppose you treat this as a problem of fitting x as a function of y? First, plot your data.
X1=[-23.24 -16.79 -14.4 -10.39 -7.494 -4.984 -2.75 0 0 0 23.24 16.79 14.4 10.39 7.494 4.984 2.75];
Y1=[-4 -3 -2.6 -2 -1.6 -1.3 -1.1 -1 0 1 4 3 2.6 2 1.6 1.3 1.1];
On top of that, I can see the curve should be visually roughly asymptotic to the line X = 6*Y. I've plotted that line in blue. And that should tell you the rough value of a, as 6.
Now, consider the model as:
if abs(y) < 1, x == 0
if abs(y) >= 1, x = a*sqrt(y^2 - 1)
this has the desired behavior. When y == +/-1, we can see that x will be zero. And it is still asymptotially linear for large y. Finally, see that I dropped b completely. Assume b==1.
Can I fit that model? Yes. I'll use the curvefitting toolbox here, the I'll do it with lsqcurvefit, since you want to use that tool.
As well, I'll show you how to write that piecewise thing using only one line, using max.
mdl = fittype('a*sign(Y1).*sqrt(max(abs(Y1),1).^2 - 1)','indep','Y1')
You should see that when abs(Y1) is 1 or larger, this function will have the shape a*sqrt(Y1^2 - 1). When abs(Y1) is 1 or less, the max function maps that to 1. But then it subtracts 1, to get 0. So this function now has the desired shape.
The sign function is NOT the trigonometric sine function.
estmdl = fit(Y1',X1',mdl,'start',1)
So your model is now
if abs(y) < 1, x == 0
if abs(y) >= 1, x = sign(y)*5.985*sqrt(y^2 - 1)
As you should see, it does a remarkably good job now. lsqcurvefit will do as easily. (Do I need to show you how to do it with lsqcurvefit too?) But don't bother trying to estimate b. You cannot do so. b == 1.