@wdarrel Applying the Chain rule (not Maple syntax):

diff(u,t) = diff(w,tau)*diff(tau,t) + diff(w,y)*diff(y,t) = diff(w,tau)*1/2 + diff(w,y)*(-1/2)

diff(u,x) = diff(w,tau)*diff(tau,x) + diff(w,y)*diff(y,x) = diff(w,tau)*1/2 + diff(w,y)*1/2

Adding those you see the cancellation of the diff(w,y) terms.

In Maple and first in general:

restart;
#The variable change:
V:={y=f(x,t),tau=g(x,t)};
eq0:=u(x,t)=w(y,tau);
eq:=subs(V,eq0);
#Differentiating:
eqt:=diff(eq,t);
eqx:=diff(eq,x);
#Finding the variable change in our case, i.e. finding y = f(x,t) and tau=g(x,t):
S:=solve({ x=tau+y,t=tau-y },{tau,y});
f:=(x,t)->(x+t)/2;
g:=(x,t)->(x-t)/2;
eqx;
eqt;
%%+%;
subs({x=tau+y,t=tau-y},rhs(%));
convert(%,diff);

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

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