@PatrickT What did you intend the last line to do?

Try

seq( Q(b,Lambda), params = [b=1,Lambda=1] ) ;

@escorpsy Don't you just mean the projection of the orbits onto the xy-plane?

If so, you can do as indicated in an earlier response, i.e. use DEplot (not DEplot3d), but add the scene option: scene=[x,y]. Beware though that there likely will be problems because some orbits fly off to infinity thereby creating lots of data with HFloat(undefined). It helps to make the timerange shorter in both ends and make the stepsize smaller (0.001).
Added: Because of the absence of similar problems with DEplot3d it may be a better idea simply to rotate the 3d plot produced by DEplot3d. This can also be done programmatically by using the orientation option. Suppose the plot generated by DEplot3d is assigned to the variable p, then
display(p,bg,orientation=[90,0,180]); #xy-plane
display(p,bg,orientation=[90,90,180]);#xz-plane
display(p,bg,orientation=[90,90,90]);##zy-plane



Another approach is to use odeplot instead as indicated in the reference given by Patrick. I haven't looked into that, but with the same definitions as before you can try:

ptsxy:=subsop~(3=NULL,pts);
bgxy:=pointplot(ptsxy,symbol=solidcircle,symbolsize=20,color=blue):
RES:=dsolve({sys,x(0)=x00,y(0)=y00,z(0)=z00},numeric,parameters=[x00,y00,z00]);
pp:=proc(x0,y0,z0,timerange)
  RES(parameters=[x0,y0,z0]);
  plots:-odeplot(RES,[x(t),y(t)],timerange, _rest)
end proc:
# pp accepts all optional arguments to odeplot (since they are passed directly to odeplot via '_rest').
pp(.4,2,1,-0.9..0,view=[-5..5,-5..5]);
PL:=[seq(seq(seq(pp(x0+i*dx/N,y0+j*dy/N,z0+k*dz/N,-1..1,view=[x0..x0+dx,y0..y0+dy]),k=1..N-1),j=1..N-1),i=1..N-1)]:
display(bgxy,PL);

@escorpsy Don't you just mean the projection of the orbits onto the xy-plane?

If so, you can do as indicated in an earlier response, i.e. use DEplot (not DEplot3d), but add the scene option: scene=[x,y]. Beware though that there likely will be problems because some orbits fly off to infinity thereby creating lots of data with HFloat(undefined). It helps to make the timerange shorter in both ends and make the stepsize smaller (0.001).
Added: Because of the absence of similar problems with DEplot3d it may be a better idea simply to rotate the 3d plot produced by DEplot3d. This can also be done programmatically by using the orientation option. Suppose the plot generated by DEplot3d is assigned to the variable p, then
display(p,bg,orientation=[90,0,180]); #xy-plane
display(p,bg,orientation=[90,90,180]);#xz-plane
display(p,bg,orientation=[90,90,90]);##zy-plane



Another approach is to use odeplot instead as indicated in the reference given by Patrick. I haven't looked into that, but with the same definitions as before you can try:

ptsxy:=subsop~(3=NULL,pts);
bgxy:=pointplot(ptsxy,symbol=solidcircle,symbolsize=20,color=blue):
RES:=dsolve({sys,x(0)=x00,y(0)=y00,z(0)=z00},numeric,parameters=[x00,y00,z00]);
pp:=proc(x0,y0,z0,timerange)
  RES(parameters=[x0,y0,z0]);
  plots:-odeplot(RES,[x(t),y(t)],timerange, _rest)
end proc:
# pp accepts all optional arguments to odeplot (since they are passed directly to odeplot via '_rest').
pp(.4,2,1,-0.9..0,view=[-5..5,-5..5]);
PL:=[seq(seq(seq(pp(x0+i*dx/N,y0+j*dy/N,z0+k*dz/N,-1..1,view=[x0..x0+dx,y0..y0+dy]),k=1..N-1),j=1..N-1),i=1..N-1)]:
display(bgxy,PL);

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

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