@User7843 The line

evalf((eval(phi,t=0)-eval(phi,t=3600))/2/Pi);

gives the number of periods of the sin function. The result is approximately 2262 on the interval I was using.
Try removing evalf. That will show what it does.
I'm using the fact that phi is a monotone (decreasing) function of t.

@Markiyan Hirnyk Quite so, thanks.

@Markiyan Hirnyk Quite so, thanks.

Events could be used instead of the loop:

ode1:=subs(A[1] = B[1], A[2] = B[2], ode):
res := dsolve({U(0)=1000, ode1}, numeric,
  events=[seq([t=k*365,U(t)=sigma*U(t)],k=1..4),[U(t)=40..1001,U(t)=0]],range=0..365*5);
plots:-odeplot(res,[t,U(t)],0..365*5,thickness=2,refine=1,color=blue);

The event  [t=k*365,U(t)=sigma*U(t)]  (k=1..4) does the essential part done in the loop.
The event  [U(t)=40..1001,U(t)=0] makes U(t) = 0 if U(t) ever falls outside 40..1001.
Notice that this is not exactly the same as is done in the loop, where the value of U(t) is only checked at the tk's.

Events could be used instead of the loop:

ode1:=subs(A[1] = B[1], A[2] = B[2], ode):
res := dsolve({U(0)=1000, ode1}, numeric,
  events=[seq([t=k*365,U(t)=sigma*U(t)],k=1..4),[U(t)=40..1001,U(t)=0]],range=0..365*5);
plots:-odeplot(res,[t,U(t)],0..365*5,thickness=2,refine=1,color=blue);

The event  [t=k*365,U(t)=sigma*U(t)]  (k=1..4) does the essential part done in the loop.
The event  [U(t)=40..1001,U(t)=0] makes U(t) = 0 if U(t) ever falls outside 40..1001.
Notice that this is not exactly the same as is done in the loop, where the value of U(t) is only checked at the tk's.

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;

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);