@zia9206314 The values of the roots computed by NeztZero in Maple 17 and Maple 18 are identical down to the last digit. It is only the computation of the residuals that is different. This is probably because of a bug in Maple 17 that was fixed in Maple 18.

No, I can't send you Maple 18.

evaluation.mw

@ecterrab I didn't say that there were several bugs; I said that I'd seen several bug reports. All the reports may have been about the same bug.

It must've been Maple 18 where that bug was fixed. Executing the worksheet evaluation.mw (from the head of this subthread) in Maple 18 and Maple 17 shows huge differences in evaluation of a HeunB function such that the residuals from both root finders (NextZero and fsolve) make no sense in Maple 17. Since the OP is using Maple 17, I think that the Heun functions can be blamed.

@zia9206314 

I could not duplicate your results. I executed your worksheet in Maple 18, and all 14 roots found by NextZero have residuals less than 5*10^(-Digits). Try executing the worksheet again. And what version of Maple are your using?

That being said, I have seen several bug reports about numerical evaluation of Heun functions. So, if there is actually a discrepancy with the residuals, I'd blame that instead of blaming NextZero or fsolve.

By the way, guardDigits= 22 is excessive. At Digits=15, guardDigits= 4 is enough to get the first 14 roots.

@Markiyan Hirnyk 

It takes guardDigits= 9 to get the first 20 roots.

I think that the OP was just blindly copying my code from an earlier problem and didn't really want 20 roots.

@ecterrab 

I don't see how evalb makes any difference. It's the default boolean evaluator anyway, and of course it's idempotent. But perhaps it's fine if it errors out from a FAIL---it's a sign that the parameters to NextZero need to be adjusted.

If your Matrix has 8 bytes per element, that would be 8*63001^2 ~ 32 GB. Does your computer have that much free memory? Since you're only using Digits = 5, you can use 4-byte floats in the matrix by declaring

m_Gi:= Matrix(num_Strat,  num_Strat, datatye= float[4]): 

Since the columns in your matrix can be generated on the fly, you may be able to use black-box methods, i.e., the columns of the matrix are generated as needed and discarded. How do you plan on using this matrix?

@shzan I already told you how to make it constant! Use convert(..., D).

Alejandro's Answer above explains the source of the problem, not how to correct it. If you are a new Maple user, I wouldn't try to understand it.

@Mac Dude I believe that the OP's problem is that the eval(diff(...)) form is not treated as a constant with respect to integration.

Compare

int(diff(f(x), x), x= -1), t)

with

int(D(f)(-1), t)

@shzan If Ans is the answer returned by the program, do convert(Ans, D).

@Markiyan Hirnyk 

That's very odd because I wrote the code in Maple 2015 Standard, mode= worksheet. What happens for you when you run it? Are you using 1D or 2D input?

What about you, Rouben? Can you run it?

@Rouben Rostamian 

By avoiding re-integrating, I got the time for your code down to 6.3 seconds.

restart:
st:= time():
Digits:= 4:
u:= Heaviside:
f:= t-> u(t) - u(t-2):
g:= t-> t*u(t) - (t-4)*u(t-4):
plotg:= plot(g, -10..10, color= blue):
h:= unapply(Int(f(t-q)*g(q), q= -10..t), t):
ploth:= convert(op([1,1], plot(h, -10..9)), listlist):
frames:= seq(
     plots:-display([
          plot(f(t-s), s= -10..10, color= blue),
          plot(select(xy-> xy[1] <= t, ploth), color= red, thickness= 3)
     ]), t= -7..9, 0.4
):
plots:-display(
     [plots:-display([frames], insequence), plotg],
     view= [DEFAULT, `..`((min,max)(ploth[..,2]))]
);
time()-st;

@nm Sorry, my mistake, I simply forgot about the extra 0.

The existence of so-called programmer indexing for rtables (A(0) instead of A[0]) makes this a little tricky. Whereas [x,y](0) becomes [x(0),y(0)], it doesn't work for <x,y>(0). We can get around that by mapping the prefix form of the function-application operator, `?()`.

Try,

map(`?()`, D~(map2(op, 0, x(t))), [0]) =~ 0;

or

<seq(D(x)(0)=0, x= map2(op, 0, x(t)))>;

I can get a series solution for both equations if I set n=1. I haven't gotten any other n to work.

@tomleslie I intended no disparagement of your Answer. You also found genuine and important bugs in the code.

The first line of code after the procedure definition applies subs to a variable U. But U has no value other than itself. That's why Tom Leslie's plots are not very interesting.