Thanks everyone for answering the question. The problem is solved.

I guess It was the problem with my expression.

@mmcdara Thank you for your response. Linearization works when r->infinity(kappa->infinity). I don't know why kappa approaches infinity( It really shouldn't). The whole point of the curve-fit was to find an expression for J(J that is not equal to 1). That's why I thought there may be a problem with my expression. However, I checked the expression last night and nothing was wrong with it. I still don't know why kappa approaches infinity.

@acer I don't know why kappa approaches infinity( It really shouldn't). The whole point of the curve-fit was to find an expression for J(J that is not equal to 1). That's why I thought there may be a problem with my expression. However, I checked the expression last night and nothing was wrong with it. I still don't know why kappa approaches infinity.

Additionally these expressions are derived from another expression containing kappa*(J-1), and I don't think kappa->infinity and J->1 is meaningful.(In physical sense J=1 means that one can ignore kappa*(J-1) which implies incompressibility.( This can't be right because the model accounts for compressibility)).

@dharr Thank you very much. This code works well. There may be a problem with my expression, because kappa approaches to infinity.

@acer Thank you for you answer. I tried 'p__eng1'(C_10, C__01, kappa, lambda) but unfortunately It didn't work. I think fsolve would be faster. I tried fsolve (r := (C__10, C__01, kappa, lambda) -> fsolve(p__3(C__10, C__01, kappa, J, lambda), J)) and I faced this error :

Error, (in Statistics:-NonlinearFit) {C__01, C__10, kappa, lambda} are in the equation, and are not solved for

This error appears even when I use 'p__eng1'(C_10, C__01, kappa, lambda).

@Ariathm P.S. I added lambda>0, because the solve creates a piecewise solution in which it turns 0 for lambda<0. I thought maybe this caused the brackets.

@nm You are right. It seems that r := (C__10, C__01, kappa, lambda) -> evalf((solve(p__3(C__10, C__01, kappa, J, lambda), J, useassumptions = true) assuming (0 < J, J::real))) produces [ ]. I added another assumption(lambda >0). but now I get another error :

Error, (in Statistics:-NonlinearFit) complex value encountered

Although I'm assuming real solutions, the complex value is encountered in the variable "r" ( I checked this by using abs() for the solve(...), which didn't cause any error).