U is not known. The output from is(U(tk)-4<=0) is FAIL (try inserting a print(is(U(tk)-4<=0));  ).

Notice the result of

if FAIL then F else G end if;

It appears that I left out a line defining 'pts':

solve({f(x, y, z) = 0, g(x, y, z) = 0, h(x, y, z) = 0}, {x, y, z});
pts:=map(subs,[%],[x,y,z]);

Another thing: In your worksheet you use 2D input and apparently somehow a space was interpreted as a multiplication sign in the definition of bg:

bg:=pointplot3d([[-1, 2/5, 0], [25/8, 5/18, -5/8]], symbol = solidsphere, symbolsize = 20, color = blue):

Advice: Don't use 2D input. I never do.

It appears that I left out a line defining 'pts':

solve({f(x, y, z) = 0, g(x, y, z) = 0, h(x, y, z) = 0}, {x, y, z});
pts:=map(subs,[%],[x,y,z]);

Another thing: In your worksheet you use 2D input and apparently somehow a space was interpreted as a multiplication sign in the definition of bg:

bg:=pointplot3d([[-1, 2/5, 0], [25/8, 5/18, -5/8]], symbol = solidsphere, symbolsize = 20, color = blue):

Advice: Don't use 2D input. I never do.

@wdarrel Make the change of variables given earlier. Then you got an ode. Solve that. Return to the original variables and you are done.

restart;
pde := diff(u(x, t), t)+diff(u(x, t), x)+u(x, t)=Dirac(x-1);
pde2:=PDEtools:-dchange({x=tau+y,t=tau-y,u(x,t)=w(y,tau)},pde,[y,tau,w]);
#First pdsolve
res:=pdsolve(pde2);
###############
#Now suppress the dependency of y and then use dsolve:
subs(w(y,tau)=wy(tau),pde2);
dsolve(%);
#_C1 will depend on y (but not tau):
subs(_C1=_F1(y),wy(tau)=w(y,tau),%);
###############
#Now back to x and t:
sbs:=solve({x=tau+y,t=tau-y},{y,tau});
subs(res,sbs,u(x,t)=w(y,tau));
res2:=combine(expand(%));
#This result is OK, but to make it look like simpler:
eval(res2,_F1=(v->F(-2*v)*exp(-v)));
combine(expand(%));

@wdarrel Make the change of variables given earlier. Then you got an ode. Solve that. Return to the original variables and you are done.

restart;
pde := diff(u(x, t), t)+diff(u(x, t), x)+u(x, t)=Dirac(x-1);
pde2:=PDEtools:-dchange({x=tau+y,t=tau-y,u(x,t)=w(y,tau)},pde,[y,tau,w]);
#First pdsolve
res:=pdsolve(pde2);
###############
#Now suppress the dependency of y and then use dsolve:
subs(w(y,tau)=wy(tau),pde2);
dsolve(%);
#_C1 will depend on y (but not tau):
subs(_C1=_F1(y),wy(tau)=w(y,tau),%);
###############
#Now back to x and t:
sbs:=solve({x=tau+y,t=tau-y},{y,tau});
subs(res,sbs,u(x,t)=w(y,tau));
res2:=combine(expand(%));
#This result is OK, but to make it look like simpler:
eval(res2,_F1=(v->F(-2*v)*exp(-v)));
combine(expand(%));

My assertion about escorpsy's system was simply that the positivity of one of the eigenvalues of the Jacobian evaluated at the critical points implied instability of the critical points. This in itself is a meager result and doesn't say much about the behavior otherwise.
And surely, pictures prove nothing. 

My assertion about escorpsy's system was simply that the positivity of one of the eigenvalues of the Jacobian evaluated at the critical points implied instability of the critical points. This in itself is a meager result and doesn't say much about the behavior otherwise.
And surely, pictures prove nothing. 

@Markiyan Hirnyk Could you please excommunicate me too?

@wdarrel I'm pretty sure that the numerical algorithms used in solving pdes in Maple do not expect distributions (here Dirac). Whatever comes out is not to be trusted.

Another point: Your pde (this and the original one) can be solved by first making a change of variables like this:
PDEtools:-dchange({x=tau+y,t=tau-y,u(x,t)=w(y,tau)},pde,[y,tau,w]);
#the result is
diff(w(y, tau), tau)+0.591e-1*w(y, tau) = 4.472135955*10^8*Dirac(tau+y-200)+1.0434983895*10^9*Dirac(tau+y-500)
You see that it has the form
(*) diff(w(y, tau), tau)+0.591e-1*w(y, tau) = f(y,tau) 
where f is known. Ignore that f is a distribution and think of it as a function of 2 variables. For every fixed y (*) is a first order ode in tau and has a unique solution if w(y,tau) is known for some tau. Since u(x,0) is known it follows that w(y,tau) is known for tau = y. Thus no more information is needed.

I suggest therefore that you proceed to solve the pde as we solved the original one. Something like this:
restart;
pde := diff(u(x, t), t)+diff(u(x, t), x)+0.591e-1*u(x, t)=447213595.5*Dirac(x-200)+1043498389.5*Dirac(x-500);
ans:=pdsolve(pde);
eval(rhs(ans),t=0)=10^8;
solve(%,{_F1(-x)});
subs(x=x-t,%);
eval(ans,%);
res:=combine(expand(%));
pdetest(res,pde);
eval(res,t=0);
eval(res,x=0);
subs(x=0,res);
U:=collect(expand(rhs(res)),exp(-591/10000*x));
h:=coeff(U,exp(-591/10000*x));
plot3d(h*10^(-21),x=0..1000,t=0..1000,axes=boxed);
plot3d(ln(U),x=0..1000,t=0..1000,axes=boxed);
plot3d(ln(U),x=0..1000,t=0..1,axes=boxed);
plots:-animate(plot,[ln(U),x=0..1000],t=0..20);

@wdarrel I'm pretty sure that the numerical algorithms used in solving pdes in Maple do not expect distributions (here Dirac). Whatever comes out is not to be trusted.

Another point: Your pde (this and the original one) can be solved by first making a change of variables like this:
PDEtools:-dchange({x=tau+y,t=tau-y,u(x,t)=w(y,tau)},pde,[y,tau,w]);
#the result is
diff(w(y, tau), tau)+0.591e-1*w(y, tau) = 4.472135955*10^8*Dirac(tau+y-200)+1.0434983895*10^9*Dirac(tau+y-500)
You see that it has the form
(*) diff(w(y, tau), tau)+0.591e-1*w(y, tau) = f(y,tau) 
where f is known. Ignore that f is a distribution and think of it as a function of 2 variables. For every fixed y (*) is a first order ode in tau and has a unique solution if w(y,tau) is known for some tau. Since u(x,0) is known it follows that w(y,tau) is known for tau = y. Thus no more information is needed.

I suggest therefore that you proceed to solve the pde as we solved the original one. Something like this:
restart;
pde := diff(u(x, t), t)+diff(u(x, t), x)+0.591e-1*u(x, t)=447213595.5*Dirac(x-200)+1043498389.5*Dirac(x-500);
ans:=pdsolve(pde);
eval(rhs(ans),t=0)=10^8;
solve(%,{_F1(-x)});
subs(x=x-t,%);
eval(ans,%);
res:=combine(expand(%));
pdetest(res,pde);
eval(res,t=0);
eval(res,x=0);
subs(x=0,res);
U:=collect(expand(rhs(res)),exp(-591/10000*x));
h:=coeff(U,exp(-591/10000*x));
plot3d(h*10^(-21),x=0..1000,t=0..1000,axes=boxed);
plot3d(ln(U),x=0..1000,t=0..1000,axes=boxed);
plot3d(ln(U),x=0..1000,t=0..1,axes=boxed);
plots:-animate(plot,[ln(U),x=0..1000],t=0..20);

@wdarrel You ned to plot the right hand side (rhs) of res and you need plot3d.

Since exponentials are involved the terms get either very large or very small, so logarithmic plots may look more interesting. Also a plot of the coefficient of exp(-x) may be illuminating:

collect(rhs(res),exp);
plot3d(Heaviside(x)-Heaviside(x-t),x=-5..5,t=0..10,axes=boxed);
plot3d(ln(rhs(res)),x=-3..3,t=0..100,axes=boxed);
plot3d(ln(rhs(res)),x=-1..3,t=0..10,axes=boxed);
#An animation
plots:-animate(plot,[ln(rhs(res)),x=-1..3],t=0..10);

@wdarrel You ned to plot the right hand side (rhs) of res and you need plot3d.

Since exponentials are involved the terms get either very large or very small, so logarithmic plots may look more interesting. Also a plot of the coefficient of exp(-x) may be illuminating:

collect(rhs(res),exp);
plot3d(Heaviside(x)-Heaviside(x-t),x=-5..5,t=0..10,axes=boxed);
plot3d(ln(rhs(res)),x=-3..3,t=0..100,axes=boxed);
plot3d(ln(rhs(res)),x=-1..3,t=0..10,axes=boxed);
#An animation
plots:-animate(plot,[ln(rhs(res)),x=-1..3],t=0..10);

Please see the additions above.

Please see the additions above.

You are getting the cpu time. The time it takes the interface to find a way to print the result is not taken into account.