@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.

@still0fthenight Many responders here are happy to help students with their homework if they have shown some effort, and you have clearly shown some effort. Now how about uploading that code?

@mapleus 

Like I said in the Answer, I needed to adjust the precision control after you filled in the signum. I've done this by using the option relerr= 1e-10 in dsolve and setting Digits:= 7 before doing the fsolve. This assumes that knowing X to seven digits is sufficient for your purposes.

Download dsolve_piecewise.mw

@itsme I'd be even more careful with the cutting and pasting. I might occasionally cut and paste from a 2d output (aka prettyprinted) region to a 1D input region (the only input I use) in order to inspect the internal structure of something (although I'm more likely to use lprint(%) for that); but I would never assume that the result of doing that is equivalent to some previous input.