@maple fan 

If you want to make them equal to 0, use evalindets: Assuming that ex is your expression and that you want to set all derivatives with respect to x to 0, do

evalindets(ex, 'diff'(anything, identical(x)), 0);

Are you looking for a numeric or a symbolic solution to your differential equations? You have three nonlinear first-order ODEs each with a little over 100 terms. Maple will try, perhaps for a very very long time and using very much memory, to solve these symbolically. Or you can get a numeric solution in a second. Either way, I know how to correct the error.

I simplified the Boolean condition by using xor. Here it is:

`&<`:= proc(x::~Vector, y::~Vector)
local n:= op(1,x), i;
     for i to n while x[i] = y[i] do end do;
     i > n or not (i::even xor x[i] > y[i])
end proc;

Both this and the previous version are intended to be run in 1D input mode. The 2D input mode has different evaluation rules for Boolean expressions. This 1D-2D difference may have been fixed in Maple 2015; I can't check that right now.

If you need high-speed performance from this procedure, let me know. The speed can be improved at the expense of the flexibility by avoiding the Vector conversion and by setting the length to a constant.

@maple fan 

Well, what do you mean by "remove"? There's only two ways that I can think of to remove a subexpression from an algebraic expression. One is to set the subexpression to 0, and the other is to set the subexpression to 1. But there might be other ways that I'm not thinking of. Please show me an example of an expression before and after the derivatives have been removed.

@tomleslie 

In any Maple program, the sequence of random numbers will be the same unless you use the command randomize().

For 1000 trials, the probabilty of deviating from 667 by more than 32 is exactly

1 - add(binomial(1000,k)*(2/3)^k*(1/3)^(1000-k), k= 667-32..667+32);

or about 3%.

For two sets of 1000 trials, assuming properly generated random numbers, the probability of getting exactly the same number of successes is exactly

add((binomial(1000,k)*(2/3)^k*(1/3)^(1000-k))^2, k= 0..1000);

or about 2%.

The standard deviation of the underlying distribution is exactly

sqrt(1000*(2/3)*(1/3));

or about 15. Thus a deviation of 32 from 667 is slightly more than two standard deviations from the mean.

For what it's worth, I can find nothing wrong with your code. It looks like you've run into a bug in dsolve. Even if your dsolve call is simplified to

dsolve(eval(sys, `&epsilon;`= L[0]), numeric)

the same error message appears. I haven't been able to figure out a workaround.

Your command is correct. To diagnose further, I need to see this "bunch of dashed lines". Can you upload a picture?

One way to upload a picture would be to right click on the graph and select Export-> PNG. Then use the green uparrow in the toolbar of the MaplePrimes editor to upload it.

@MTata94 

1-D Math Input means that the characters that you type are exactly what occurs in the expression, and they occur from left to right at your cursor. Most (not all) of the oldtimers here at MaplePrimes use 1-D Math Input in one of the three User Interfaces (UIs) offered by Maplesoft: Standard, Classic, or Command Line. I use 1-D with Standard. To select this, go to the Tools-> Options menu. Select the Display tab. In the first pull-down "Input Display", select "Maple Notation". Go to the Interface tab. In the fourth pul-down "Default format for new worksheets", select "Worksheet". The click Apply Globally. Then open a new worksheet.

@85951971 

If you want to select the model based on the lowest residuals, it is possible to do that automatically, but not necesarily advisable. It is very tricky with nonlinear models: There are often parameter choices that give a local minimum for the residuals which are far from being a global minimum. Maple's Statistics:-NonlinearFit (and all nonlinear parameter fitting software that I am aware of) very often returns one of these local minima.

Take a close look at Markiyan's answer. The program that he used is called DirectSearch:-DataFit (the DirectSearch package must be downloaded from the Maple Applications Center). His correct piecewise interpretation of the spherical model gives a good fit. Notice the first number in the output of DataFit. That is the total residual. From all the fitmethod choices offered by DataFit, Markiyan hand selected the one with the smallest first number. It would be easy to automate that.

Markiyan (who is from Ukraine) has an excellent knowledge of English. If you write in correct English with correct spelling and punctuation and without hipster slang he will understand you.

@Rouben Rostamian  

Thanks Rouben. It's very impressive.

@Preben Alsholm 

If I simply differentiate eq2 and use a DAE method, then I get the same solution as your "Alternative solution". It would seem like the loss of information caused by differentiating the constant would make this method not work, and yet the differences between the two numerical solutions are insignificant (on the order of 10^(-8)).

I simply used

solution2:= dsolve([diff(eq1, t), eq2, CI], numeric, method= rkf45_dae);

Do you have any idea why this works?

@MDD In other words, the list L of monomials needs to be generated automatically, and needs to be in the correct order.

@acer I updated the Answer. I can find no significant difference in the times when hits are unlikely, and, yes, there are other factors to consider that my timing test does not address. I understand and appreciate your argument based on clear coding style.

In the more likely case that hits are likely, I found that the assigned method is significantly faster.

@Markiyan Hirnyk 

Yes, it is not necessary to assign the plot to a variable and then display it. Simply using the odeplot command, without assignment, is enough.

@Wolff 

The dialog box can be manually moved to reveal the whole plot. But that's a pain to do every time.

Here's a completely different idea: Put the plots in an animation, and step through them in single-frame mode, like a slideshow:

Test:= proc(j)
local x;
     plot(x^j, x= -2..2)
end proc:

plots:-display([seq(Test(k), k= 1..3)], insequence);