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

@escorpsy I have uploaded a worksheet containing illustrations and animations. No text, but hopefully somewhat understandable.

MaplePrimes12-11-2.mw

@escorpsy I have uploaded a worksheet containing illustrations and animations. No text, but hopefully somewhat understandable.

MaplePrimes12-11-2.mw

@Markiyan Hirnyk I gave you Brannan and Boyce as a general reference since it wasn't clear to me how much you knew about these things. But you are right, B&B doesn't cover the present situation.

You may want to take a look at

http://www.scholarpedia.org/article/Chetaev_function

After changing variables the present system can be written in vector form as

dx/dt = Dx+|x|e(x), where e(x) ->0 as x->0 and D is a diagonal matrix having the eigenvalues in the diagonal, say in the order -1, 0, 6.

As a Chetaev function you can take V(x1,x2,x3) = -x1^2 + 6*x3^2.

I leave it as an exercise for you to work out the details.