Sadly, while you want some magical form for each of these curves, that is not so easy to find. Essentially, it looks like each of those curves is well modeled by a simple low order polynomial model, coupled with a maximum, a saturation point. Nothing more sophisticated is necessry to do the above. And there is not sufficient information content in the data to build somethign more sophisticateed anyway.
y1=[2.1 3.045 4.2 4.515 6.09 6.72 6.8145 6.825 6.825 6.825 6.72 6.825 6.825 6.93 6.825 6.825];
y2=[1.89 2.31 2.625 2.94 3.675 3.99 4.41 5.04 6.09 6.51 8.190 9.345 10.5 10.92 10.92 10.92 10.92 10.71 10.92 10.5 10.92 10.92];
y3=[2.415 2.625 3.045 3.675 3.99 4.2 4.83 5.565 5.88 6.51 6.825 7.875 8.61 9.975 10.5 11.34 11.55 13.965 14.9625 15.12 15.12 15.33 15.12 14.805 15.12 14.91 15.12 15.12 15.225 15.12];
y4=[2.94 3.15 3.255 3.57 3.885 4.305 4.41 5.04 5.25 5.67 6.09 6.3 6.615 7.77 7.98 8.631 9.24 9.975 10.605 11.34 12.075 13.335 14.175 14.7 16.17 17.325 18.585 19.215 19.2465 19.215 19.215 19.215 19.215 19.425 19.215 19.215 18.9 19.215 19.215 18.69 19.215 19.215 19.2150];
If we look at the above curve, we can see a straight line, coupled with a saturation level. We might do that iusing a simple model, such as...
Linmdl = fittype('a + b*min(x,xsat)','indep','x')
Linmdl =
General model:
Linmdl(a,b,xsat,x) = a + b*min(x,xsat)
So xsat is the point at which the curve saturates. All that really matters is I choose at least a decent estimate of the saturation point location.
fittedmdl1 = fit(x1',y1',Linmdl,'start',[1 1 2])
fittedmdl1 =
General model:
fittedmdl1(x) = a + b*min(x,xsat)
Coefficients (with 95% confidence bounds):
a = -0.67 (-1.116, -0.2243)
b = 3.72 (3.41, 4.03)
xsat = 2.015 (1.945, 2.084)
There simply is nothing more in that data. No sophisticated model using tan or exp. Your data justifies nothing more. The y value at the saturation point is simple to find too.
ysat = fittedmdl1(fittedmdl1.xsat)
Similarly, we can model each of the other curves the same way. There, I could use a quadratic polynomial, since there is some amount of curvature in the curve below the break point.
quadmdl = fittype('a + b*min(x,xsat) + c*min(x,xsat).^2','indep','x')
quadmdl =
General model:
quadmdl(a,b,c,xsat,x) = a + b*min(x,xsat) + c*min(x,xsat).^2
fittedmdl2 = fit(x2',y2',quadmdl,'start',[1 1 1 4])
fittedmdl2 =
General model:
fittedmdl2(x) = a + b*min(x,xsat) + c*min(x,xsat).^2
Coefficients (with 95% confidence bounds):
a = 2.271 (1.557, 2.984)
b = -0.776 (-1.478, -0.07389)
c = 0.7863 (0.6328, 0.9397)
xsat = 3.833 (3.763, 3.904)
ysat = fittedmdl2(fittedmdl2.xsat)
Again, we see a very nice fit, entirely reasonable. I'll do a fit for the third curve, just to show how easy it is. Again, I'll use the quadmdl form.
fittedmdl3 = fit(x3',y3',quadmdl,'start',[1 1 1 6])
fittedmdl3 =
General model:
fittedmdl3(x) = a + b*min(x,xsat) + c*min(x,xsat).^2
Coefficients (with 95% confidence bounds):
a = 2.205 (1.336, 3.074)
b = -0.2222 (-0.7667, 0.3223)
c = 0.4144 (0.3378, 0.4911)
xsat = 5.853 (5.763, 5.942)
You can fit the last curve. A quadratic model will surely suffice there too for the lower region. Again, if you think some complicated model with exp and tan functions is necessary, you are simply wrong. The data does not justify it, nor could you estimate such a model. And since you provide no physical reason why some specific form should exist, we cannot possibly guess what it might be.
DON'T OVERFIT YOUR DATA! Apply basic common sense to what you have, and don't look for something, more complicated. While that more complex model might exist in the shadows, you won't be able to know what it is. And if you did try to find something complex, you would get complete garbage as a result of that exploration on such limited information.