Sivan,

1. I think that you might be missing a comma. Are your equations:

   {P1'[t] == -(2 y + x) P1[t] + (x - y) P2[t] + y P2'[t] == -(2 z + x) P2[t] + (x - z) P1[t] + z}
   

or

{P1'[t] == -(2 y + x) P1[t] + (x - y) P2[t] + y, P2'[t] == -(2 z + x) P2[t] + (x - z) P1[t] + z}

I believe you probably meant the later but you coded the former (both are valid although the first, the one you are using, is an unusual way to express the system of equations).

1. you are doing something computationally evil. You are recomputing the solution to a differential equation over and over for different data. The differential equation should be solved once and reused.

For example:

nsol = DSolveValue[{P1'[t] == -(2 y + x) P1[t] + (x - y) P2[t] + y, 
     P2'[t] == -(2 z + x) P2[t] + (x - z) P1[t] + z, P1[0] == 1, 
     P2[0] == 0}, P1[t] - P2[t], t] // FullSimplify;

This will give you the response you use (P1-P2) as a function of x,y,and z.

1. why are you fitting an exponential to your solution? The solution to a second order linear differential equation is the summation of two exponential terms. These two terms will never change and you can get them analytically as a function of x,y,and z. What is the meaning of your blended exponential answer? (Of course trying to fit more than one exponential is numerically unstable so I would not try that). If you explain your goal, I can make a better suggestion.

Regards,

Neil