I have a system of six, linear ordinary differential equations. There are six parameters to be estimated in the system. It's nothing fancy. I have six variables for which I have the data from an experiment. I want to find the best fit to the data by estimating the values of the six parameters. To make matters just a bit more difficult, the ode integration has to start at time t=18 hours because the time 0 data point is not considered as relevant to the fitting process. So here is the matrix for the odes. The variable names are respectively u,v,w,x,y,z (t).

MatrixForm[{{-b - d0, 0, 0, 0, 0, 0}, {0, -bp - d0p, 0, 0, 2 a3 p, 
   0}, {b, bp, -a1 - d1, 0, 2 a3 (1 - p), 0}, {0, 0, a1, -a2 - d2, 0, 
   0}, {0, 0, 0, a2, -a3 - d3, 0}, {d0, d0p, d1, d2, d3, -dd}}]

The data can be formatted as needed but currently is formatted as follows:

data = {
   {18.00, 3483.50, 0.00, 392.88, 26.19, 0.00, 550.03},
          {24.00, 3686.30, 7.10, 312.52, 46.17, 7.10, 990.83},
          {42.00, 2409.23, 31.82, 1118.24, 877.32, 718.22, 2309.22},
          {48.00, 1660.65, 44.73, 1045.60, 1157.42, 944.95, 1990.54},
          {66.00, 293.58, 55.05, 2642.26, 1100.94, 1146.81, 4963.41},
          {72.00, 384.18, 72.03, 3457.62, 1440.67, 1500.70, 6495.04}
   };

where the first value in each group is the time (hours) and the next six values are the measured experimental values of u,v,w,x,y and z(t). The relevant ode system is given as

pfunfull = ParametricNDSolveValue[
   {u'[t] == -(b + d0) u[t],
    v'[t] == 2 p a3 y[t] - (bp + d0p) v[t],
    w'[t] == b u[t] + bp v[t] - (a1 + d1) w[t] + 2 (1 - p) a3 y[t],
    x'[t] == a1 w[t] - (a2 + d2) x[t],
    y'[t] == a2 x[t] - (a3 + d3) y[t],
    z'[t] == 
     d0 u[t] + d0p v[t] + d1 w[t] + d2 x[t] + d3 y[t] - dd z[t],
    u[tstart] == 
\!\(\*SubscriptBox[\(data\), \(\(\[LeftDoubleBracket]\)\(1, 
      2\)\(\[RightDoubleBracket]\)\)]\) ,
    v[tstart] == 
\!\(\*SubscriptBox[\(data\), \(\(\[LeftDoubleBracket]\)\(1, 
      3\)\(\[RightDoubleBracket]\)\)]\),
    w[tstart] == 
\!\(\*SubscriptBox[\(data\), \(\(\[LeftDoubleBracket]\)\(1, 
      4\)\(\[RightDoubleBracket]\)\)]\),
    x[tstart] == 
\!\(\*SubscriptBox[\(data\), \(\(\[LeftDoubleBracket]\)\(1, 
      5\)\(\[RightDoubleBracket]\)\)]\),
    y[tstart] == 
\!\(\*SubscriptBox[\(data\), \(\(\[LeftDoubleBracket]\)\(1, 
      6\)\(\[RightDoubleBracket]\)\)]\),
    z[tstart] == 
\!\(\*SubscriptBox[\(data\), \(\(\[LeftDoubleBracket]\)\(1, 
      7\)\(\[RightDoubleBracket]\)\)]\)}, {u, v, w, x, y, z}, {t, 
    tstart, tend},
   {b, d0, p, a3, bp, d0p, a1, d1, a2, d2, d3, dd}
   ];

I have read the various community work notes and gone online and even talked with Mathematica tech support. Tech support says you cannot use pfunc because I have 6 variables and parameters that I am simultaneously trying to fit and to estimate. When I asked them how to do it in Mathematica, I got the following very helpful answer: "This will be up to you."

I am not a Mathematica wizard or even a geek. I could use any help that anyone has to offer. And if you want to communicate directly, you can reach me at tmwitten@vcu.edu ... Thank you everyone. I attach the notebook with the program.

Attachments: