How do I fit experimental data by solving 3 coupled ODEs with adjustable parameters?

I want to fit experimental data by solving 3 coupled differential equations.The ODEs are rate equations modeling the dynamics of hot carriers in semiconductors after excitation. The function use fot fitting data is x2 and the parameters are ntrmax, trec and tc. I was inspired by comments from Peng Bai when writting the following code but It does not work and I can't figure out why yet.

P[t_]= 4.9560 * 10^18 Exp[-1732.87 t^2];

Nphoton = 2.1102 * 10^17;
tr = 0.2;

data = {{{0., 0., -0.000367822, 0., 0.000435918}, {0.1334, 0.01008,0.0000753379, -0.000200166, -0.000520431},{0.267,0.02302,-0.000766664,-0.000244647,-0.000209062},{0.4004,0.03347, 0.00207, 0.00133,-0.0000533769},{0.5338,0.06225,0.00772, 0.00654, 0.00186},{0.6672, 0.07971, 0.01499, 0.01052,0.00362}, {0.8006, 0.07766, 0.0196, 0.01717, 0.00522},{0.934,0.07577, 0.03668, 0.02662, 0.00975}, {1.0674, 0.07153, 0.03945,0.02604,0.01389},{1.2008,0.0715,0.03943,0.02398,0.01189}, {1.3344, 0.06899,0.03748,0.02269,0.01147},{1.4678,0.06541, 0.03571, 0.02017, 0.01113},{1.6012,0.06087,0.03342, 0.01928, 0.01104},{1.7346, 0.06103, 0.03238, 0.01886, 0.01029}, {1.868, 0.06061, 0.02904, 0.0179, 0.00922},{2.0014,0.05896, 0.02939, 0.01806, 0.00771}, {2.1348, 0.05629, 0.0285,0.01657, 0.00849},{2.2684,0.05387, 0.02817,0.01579,0.00782}, {2.4018, 0.05418, 0.02622, 0.0157, 0.00693}, {2.5352,0.05342,0.02558,0.01506,0.00633},{2.6686, 0.05169, 0.0235,0.01417, 0.00626}, {2.802, 0.04953, 0.02494,0.01348,0.00568},{2.9354,0.04977,0.02392, 0.01312, 0.00528},{3.0688, 0.04788, 0.02347, 0.01241, 0.00573},{3.2022, 0.04753, 0.02294,0.01152,0.00609},{3.3358, 0.04713, 0.02148, 0.01212,0.00488},{3.4692, 0.04553, 0.0213, 0.01136, 0.00504},{3.6026,0.04461, 0.02146, 0.01079,0.00504},{3.736, 0.04348, 0.02108,0.01025, 0.00455},{3.8694, 0.04279,0.01969, 0.01021,0.00419},{4.0028, 0.04337, 0.01864, 0.01021,0.00522},{4.1362,0.04321, 0.01913, 0.01003, 0.00466},{4.2696,0.04168, 0.01842,0.00914, 0.00448},{4.4032, 0.0395, 0.01758, 0.00952,0.00442},{4.5366,0.03814, 0.01689, 0.00843, 0.00399},{4.67,0.0403, 0.01647, 0.00854,0.00391}, {4.8034, 0.03894, 0.01654,0.0083, 0.00364},{4.9368,0.03865, 0.01654, 0.0091,0.00397},{5.0702, 0.03732, 0.0168, 0.00896, 0.00233},{5.2036, 0.03672, 0.01696, 0.00812, 0.00364},{5.337, 0.0359, 0.01592, 0.00825, 0.00344},{5.4706, 0.03541, 0.01621, 0.00823, 0.00322},{5.604, 0.03381, 0.0157, 0.0083, 0.00397},{5.7374, 0.03472, 0.01514, 0.00747, 0.00275},{5.8708, 0.03456, 0.01525,0.00754, 0.00368},{6.0042, 0.03414, 0.01454,0.00705,0.00342},{6.1376, 0.03198, 0.01472, 0.00656, 0.00266},{6.271,0.03209, 0.01421, 0.00705, 0.00326},{6.4044, 0.03245, 0.01355,0.0079,0.00315},{6.538, 0.0312, 0.01386, 0.00696, 0.00295},{6.6714, 0.03103, 0.01377, 0.00623, 0.00213},{6.8048, 0.03118, 0.01304, 0.00643, 0.00264},{6.9382, 0.03016, 0.01135, 0.00718, 0.00262},{7.0716, 0.02962,0.01189, 0.00718, 0.00159},{7.205, 0.02998, 0.01189, 0.00585, 0.00244},{7.3384, 0.02956, 0.01171, 0.00607, 0.00226},{7.472, 0.02885, 0.01246,0.00621, 0.00233},{7.6054, 0.02869, 0.01235, 0.00616,0.00239}, {7.7388, 0.02882, 0.01162, 0.00614, 0.00235},{7.8722, 0.02831,0.01233,0.00583, 0.00255}, {8.0056, 0.02664, 0.01111,0.00634,0.0027}, {8.139, 0.02747, 0.01038, 0.00518,0.00295}, {8.2724, 0.0262, 0.0106, 0.00589,0.00264}, {8.4058,0.02662, 0.01009, 0.00547, 0.00226}, {8.5394,0.0256, 0.00947, 0.00527, 0.00206}, {8.6728, 0.02524, 0.0102,0.00467,0.00206}, {8.8062, 0.02415, 0.00971, 0.00458, 0.00197}, {8.9396,0.02302, 0.00945, 0.00434, 0.0019}, {9.073, 0.02382, 0.00969, 0.005, 0.0021}, {9.2064, 0.02346, 0.00989, 0.00527,0.00175}, {9.3398, 0.02289, 0.00967, 0.00483, 0.00179}, {9.4732,0.02226, 0.01007, 0.00527, 0.00177}, {9.6068, 0.02291, 0.00925, 0.00451, 0.00133}, {9.7402, 0.02264, 0.00985, 0.0052,0.00204}, {9.8722, 0.02362, 0.00909, 0.00405, 0.00104}, {10.007,0.0226, 0.00954, 0.00451, 0.0013}, {10.1404, 0.02137, 0.00907,0.0044, 0.00199}, {10.2738, 0.02186, 0.00958, 0.00451, 0.00195}}}

ti = data[[1, All, 1]];
ci = Nphoton*data[[1, All, 2 ;; 5]];

pfun = ParametricNDSolveValue[{x1'[t] == P[t] - x1[t]/tr,x2'[t] == x1[t]/tr - x2[t] (1 - (x3[t]/ntrmax))/tp - x2[t]/trec, x3'[t] == x2[t] (1 - (x3[t]/ntrmax))/tp - x3[t]/trec, x1[0] == P[0], x2[0] == 0, x3[0] == 0}, {x1, x2, x3}, {t, 0, 20}, {ntrmax, trec, tp}];(*three dependent variables*)

f[ntrmax_?NumericQ, trec_?NumericQ, tp_?NumericQ] := Sum[Total[(ci[[All, i]] - Map[pfun[ntrmax, trec, tp][[i]], ti])^2], {i, 1, 4}] // Quiet;

fit = NMinimize[f[ntrmax, trec, tp], {ntrmax, trec, tp}];

params = fit // Last

• In your definition of f your pfun only has 3 interpolating functions and you are taking parts {i, 1, 4} which gives a string of "Part 4 of {<<1>>} does not exist" errors. Changing that 4 to 3 seems to get rid of those errors, but you have another error which may be that you aren't correctly getting a numeric result from pfun. Can you track down and correct those? You also have extremely large and extremely small coefficients and only a single digit of precision for tr, all of which may make for very uncertain results. Without trying to do the fit, can you get good results from f?
– Bill
Commented Aug 18, 2018 at 19:32
• data has five columns. What do these columns correspond to among the parameters t, x1, x2, x3? Commented Aug 19, 2018 at 2:32
• @BilI I changed the number of interpolating functions to 3 in pfun tried to simply the problem by normalizing the data ci = data[[1,All, 2;;5]] and pfun = ParametricNDSolveValue[{x1'[t] == P[t]/P[0] - x1[t]/tr,x2'[t] == x1[t]/tr - x2[t] (1 - (x3[t]/ntrmax))/tp - x2[t]/trec, x3'[t] == x2[t] (1 - (x3[t]/ntrmax))/tp - x3[t]/trec, x1[0] == 1, x2[0] == 0, x3[0] == 0}, {x1, x2, x3}, {t, 0, 20}, {ntrmax, trec, tp}]; Computed parameters do not converge to experimental data. Commented Aug 19, 2018 at 3:02
• The first column (left) in data refers to independent variables and three columns correspond to dependent variables. Commented Aug 19, 2018 at 3:07