You have input data that is not described at all. The names used are not at all mnemonic (presumably Xtime is time in minutes and Yconcentration is concentration is some unspecified units. You should have this described in your sheet, if the sheet is to be useful for anything more than a one Xtime use, to be then discarded and never looked at again. While Mathcads ODESolvers make using units a pain, you should either use units in Mathcad, or at least document the units for all variables. To further confuse things, here you label the Xtime as time and the concetration as Y. But later you use t or Xtime for time, and use x for the mass quantity. Consistent and mnemonic nomenclature helps make a sheet clearer, and reduces the chances of error.
You seem to be looking to see the extent to which the late samples affect the result. I'm not sure of the usefulness of that, a plot of the residuals (which should always be done) will show if there are anomalies with the late data. But still .....
You should be much more careful with your units, and the meaning of the variables. More definitions (verbal) are useful. Here you describe c10 and c20 as concentrations, but you eventually use them as initial values for the x1 and x2 variables, which are mass quantities. That is not consistent. Since the values are zero it doesn't affect the evential answers, but it is still logically incorrect.
Again, more description and documentation needed. This is the external infusion rate (by inference). f is not very mnemonic for this function.
f(time) is the infusion rate and could be a konstant dose over the infusion time (= the duration of the injection, t) and after the t it could be a constant rate, when a pump is used. Otherwise it is 0. In our case we do not use a pump so it is r=0.
You are trying to do too much in this structure. You do not actually need the residuals, the predicted values are quite sufficient. But you do want to use the solution both for the parameter estimation and also for plotting. You should not have a duplicate of the ODESolve block for the two purposes. That's what functions are for, to reuse the same definition for multiple calculations.
There are still some issues with this. You have not really described what the variables are. x1 and x2 are the mass quantities in the two compartments. That should be stated explicitly. Especially as you are also working with concentrations (you have expressed the initial values as concentrations, and your input data is a concentration). You could set up the equation to be in terms of concentrations rather than mass quantities. Then you would have only one rate constant, but the two volumes would be involved in the equations. As it is the volumes are implicit in the rate constants and do not explicitly appear in the ODE. The volume is not needed as a parameter.
Eliminationrate k01*Central Volume = CLEARANCE [ml/min]
t21[min] = Half life of tranfer Vol2 to Vol1 = ln(2)/k21
k21 is the constant for transfer from V1 to V2. It is not directly related to elimination. So t21 is the half life for equilibriation between the two compartments, not the half life of elimination.
Having done this once, it is unecessary to do it again. Besides which, as written, it is wrong. Sol is described with arguments of a, b, c and d, but these variables do not appear in the definition and so are meaningless. The values of the variables used are those defined above.
The Maple code is a mess. It is full of variables that I have no idea what they mean or how they were calculated. There is nothing to relate the various expressions there to the original ODE. First, I will ignore all that and just look to duplicate the Maple solution in Mathcad.
The Maple stuff seems to be set up in terms of vol1, vol2, t21, and clearance. While a bit wordy, I will stick to that.
A variant on this is to combine all the calculations into a single function. This allows using local variables rather than function definitions for the various intermediate variables. The advantage is much smaller, and more readable, expressions -- close to the Maple expressions. The disadvantage is that there is no easy way to evaluate and see the intermediate variables, making it difficult to validate the construct.
If you are going to analyze partial data, the place to set it up is here, rather than duplicate code later:
Mixing Matrix: In the rows there is the same substance and in the columns there is the mixingrelationship.
Average artificial Protocol based on M1 and M2, on these models the Vol2 is changed.
(In the case of 2 changes: M1 has the Vol2=32.3720 l and M2 has the Vol2=10.2262 l)
To use this for different analyses, you need to porivde the endpoint and also parameterize the parameter values (rather than use the current worksheet values).
Central differences are used for changing the parameters:
This looks for the first Parameter (vol1) as follows:
Parameterchange for vol1 left from the central point
Parameterchange for vol1 right from the central point
Since the residuals are often of interest in themselves, I prefer an approach that calculates the residuals and then uses them directly for the standard error calculations, rather than doing an explicit sum of squares.
For the Monte Carlo Simulation we need a interferenced curve with a normal distributed random number vector:
The discussion was the question, if you must interference the random numbers with the idiale curve (that is the Estelberger assumption) or with the observed Yconcentration values. (that is the wrong assumption)
Guess values can be provided by worksheet assignments, or as arguments of the function. There is no need to do both. Since runs are for quite similar values, there is no need to parameterize the guess values, fixed values will do fine.
SParaA is exactly MINI, without the extraneous parameters. Since the extraneous parameters were removed, it is exactly the same.
It's taking entirely too long to calculate SysParaArti many, many times. Once for each parameter value is enough!
RESULTS of system parameters with standard deviations from Fisher's Info Matrix (central differences):
Note that different data sets (truncations) have different parameters and so the model and the residuals need to be calculated with the current values.
Please notice: The assumption that the original Yconcentration should be nearly the same with the YMmean is not correct. You can see this, that the dlearance of the original concentration increases with time, and that is not with the YMmean concentration. So the mixing of the concentrations are not the same with the reality, as in reality the small molecules are feed into the concentration continously.
Column 0 is for the Systemparameters like in My3DConv above.
Column 1 is for the standard derviations done by Fisher`s sensivity analyses
Column 2 is the Fisher Info matrix.
Column 3 is the corresponding t21 value.
My preference for an analytic solution is one that starts with the actual ODE and is traceable back to that.
A simple general solution requires the exponential of a matrix. Easy enough to do with the eigenvector decomposition of a matrix. But what we really need is 
for a given A and many values of t. Eigenvector calculations are somewhat time-consuming. So we design a function that can, given A, create the function 
as a function of t.
Since we want to use the symbolic processor to extract information for the ODEs, we need to symbolically undefine variables that we wish to have symbolically.
We write the ODEs in the normal style, without the explict functional dependencies. We store them in a variable for further use.
We get complaints from the numeric processor about undefined variables. Those errors don't matter, the symbolic processor deals with undefined variables just fine. But the numeric processor gets into a huff when it is ignored. Just ignore it.
Now we do a bit of general calculus. Consider the differential equation 
where A and B are constants. The solution will be of the form 
where C and D are constants. We can find D by differentiating the equation for y and equating it to the value from the ODE, and then equating terms. We get 
. We can solve that last equality for D (the y's cancel out) and get 
. We can find C from the inital (or boundary) value. If we have that 
then we can substitute and solve to get 
. Using these two solutions we can write a function to get C and D:
We can get the matrix A from the derivative expression by taking the Jacobian. A is a function of the three k values.
We can now write a function to generate the analytic solution. But we note that B is only piecewise constant, having one value from 0 to t and a different value beyond that. So we need to generate the solution as a piecewise function.
A modified model assuming a non-zero r applied to compartment 2, and an initial value corresponding to the equilibrium.
