@fatemeh1090 

First, assign the dsolve results to a variable, as

sol:= dsolve(...);

Then, to plot p(t), use

plots:-odeplot(sol, [t, p(t)], t= 0..1);

and do likewise for the other functions. Ordinarily, I'd put the functions together in a single plot, but it's not possible in this case due to extreme differences in the scaling of the vertical axis.

@HS There would need to be some fairly special relationship among the coefficients for the gcd to not be 1. 

@eggs1996 If you have an improved procedure that you'd like further help with, let me know. 

Are you expecting a result other than 1, such as x^3 - e where e is some complicated expression that depends on a, b, and c?

@syhue Yes, my main procedure works for any finite number of dimensions. The extra plotting procedures are for 2-D data.

@666 jvbasha If I execute your worsheet top-down, I get "Error, (in evalf/int) invalid arguments". This suggests that the output that you show is not the direct result of the code that you show. That could possibly be the result of executing things out of order.

@mmcdara I included the column for your "reduced uncertainty", but I just called it "z-score" because I think that's what it's more widely called (at least in English). I also kept the p-value. I changed the **** flag to a variable-length number of stars (0 - 5) depending on how low p is. I added an option to test either the left or right side (with the default being the right), which is important for the nonsymmetric distributions. You can see all this in my new Answer.

@epostma Hi Erik,

Please see my Answer below where I've tested every stock continuous distribution for this sampling irregularity. The ones that have it are ChiSquare, Error, FRatio, InverseGaussian, LogNormal, and Normal. The Exponential distribution does not show this anomaly, so I guess that the problem is in the "fallback algorithm for the tail", which isn't needed for the Exponential.

Let me know if you have any questions about my testing method.

@mmcdara What you call "reduced uncertainty" is not normalized to N^[*1], making it difficult to put into context. I think that it should be converted to p-values as commonly used in hypothesis testing of proportions. I've done this, making the code:

SampleCheck:= proc(X, N::posint)
uses St= Statistics;
local
    p:= 10.^~[$-6..-1],
    S:= St:-Sample(X, N),
    Obs:= rhs~(St:-TallyInto(S, `..`~(min(S)-1, St:-Quantile~(X, p, 'numeric')))), 
    Exc := N*~p,
    #Compute two-tailed p-values using Normal approx to binomial:
    U:= 1 -~ abs~(1 -~ 2*St:-CDF~(Normal~(Exc, sqrt~(Exc*~(1-~p))), Obs))
;
    printf("  prob      exact    observed           p-value \n");
    printf~(
        "%1.1e  %8d %11d           %0.4e %s\n", 
        p, round~(Exc), Obs, U, `if`~(U <~ 0.05, "****", "")
    );
    return
end proc
:

SampleCheck(Statistics:-RandomVariable(Normal(0,1)), 9*10^7);
  prob      exact    observed           p-value 
1.0e-06        90         353           0.0000e+000 ****
1.0e-05       900        2016           0.0000e+000 ****
1.0e-04      9000       10435           0.0000e+000 ****
1.0e-03     90000       90131           6.6219e-001 
1.0e-02    900000      899403           5.2708e-001 
1.0e-01   9000000     8996536           2.2356e-001 

[*1] Edit: Upon taking a closer look, I realized that the problem was that you made a mistake computing the standard deviation of the binomal. It should be sqrt~(Exc*~(1-~p)), whereas you had sqrt~(Exc). With this change, your "reduced uncertainty" makes sense. Still, I prefer to further process it to a p-value.

@alexk213 To do your symbolic sum, I found that it worked better to change piecewise to `if` and Integrate to int:

phi:= (j,x)-> `if`( #works better than piecewise in this case
    j=0, exp(x),
    int(exp((1-theta)*x)*theta^(j-1),theta=0..1)/(j-1)! #Note: int
):
phi(0,x)/x^j - sum(phi(i,0)/x^(j-i), i= 0..j-1);

The result of this is

exp(x)/x^j - exp(x)/x^j*GAMMA(j,x)/GAMMA(j)

Clearly, this is equivalent to 

exp(x)/x^j*(1 - GAMMA(j,x)/GAMMA(j))

although it's tricky to get Maple to put it in this exact form.

@sand15 In my opinion, Maple 2019 is the best version since objects were introduced (Maple 16?).

There was a major bug in the 2019.2 update (restart not reading initialization file). They very recently announced that it's been fixed. 

@alexk213 Sometimes to recompute a result at a higher precision, it is necessary to "clear" the lower-precision results from the remember tables. The easiest way is restart. More subtle clearing can be done with the forget command.

I can't confirm that this is the actual source of your anomaly without seeing your code and knowing your Maple version. Like the other respondents, I can't reproduce the anomaly.

@acer I suspect that the OP was simply unaware of the formula that the transpose of a product equals the reversed-order product of the transposes of the factors.

@Kitonum Your command works fine for this unusual case, and, indeed, it makes sense intuitively. However, it should be pointed out for unwary readers that if the input had been an actual series structure (as returned by series or taylor), which is the usual case, the remove command wouldn't work.

@acer I think that the OP has a good point in this case because 1) convert(..., polynom) is a confusing name and 2) this specific case is not a series structure.

A better convert command for converting series would be convert(..., `+`), but it doesn't work.