That error message manifests a "misunderstanding" by Maple but it appears to be harmless.  Introduce a third parameter, k, and assign it an arbitrary value and Maple produces the correct result:

restart;
f:=x->piecewise(x>=-1/2 and x<=1/2,1);
eq:={diff(x(t),t)-sum(f(t-k*T),k=0..K)*x(t)=0,x(0)=1};
sol:=dsolve(eq, numeric, parameters=[T,K,k]);
sol(parameters=[1,3,-1]);  # the third argument is an arbitrary number
plots:-odeplot(sol, [t, x(t)], t=0..6,
    color=red, thickness=3, view=0..40);

Or with a more interesting choice of parameters:

sol(parameters=[2,2,-1]);  # the third argument is an arbitrary number
plots:-odeplot(sol, [t, x(t)], t=0..8,
    color=red, thickness=3, view=0..13);

Download mw.mw

Write the summation as

sum( irem(s+r,2) * whatever, s=1 .. r-1);

because irem(s+r,2) is zero when s+r is even, and one otherwise.

Asking for a phase portrait for your equation is not exactly meaningful.  The responses that you have received present graphs of solutions y(x) versus x, which is not the usual meaning of a phase portrait.  A phase portrait is a plot in the phase space.

The concept of the phase space, however, pertains to autonomous differential equations.  Your equation is non-autonomous, so the graphs that you have been presented with are remotely reminiscent of phase portraits but you should be aware that they are not phase portraits.

That said, here is what I would do if I were to plot a set of graphs of y(x) versus x (not a phase portrait).

>  restart; >  with(DEtools): >  de := diff(y(x),x)=x^2-y(x)^2; >  a := 2.0;  b := 3.0; >  ic := seq(y(s)=+a, s=-a..a, 0.3),       seq(y(s)=-a, s=-a..a, 0.3),       seq(y(-a)=s, s=-a..a, 0.3),       seq(y(+a)=s, s=-a..a, 0.3): >  DEplot(de, y(x), x=-b..b, [ic], y=-b..b,     linecolor=black, thickness=1, arrows=none);

Download mw.mw

Download mw.mw

Define φ and ψψ as you already have.

Then:

r := 2^k * M;
phi_vec := Vector[column](r, i-> phi[i]);
`&psi;&psi;_vec` := Vector[row](r, j-> `&psi;&psi;`[j]);
P := phi_vec . `&psi;&psi;_vec`;

Aside: The symbol `&psi;&psi;` in your calculations is too cumbersome and unnecessary. Why not use a single-letter symbol such as an uppercase Ψ insteaad?

Here is a worksheet that shows the general method of assembling a matrix from smaller matrices.

I must alert you, however, to your use of M^(−1) *~ K in your worksheet. That certainly IS NOT the same as the matrix product M^(−1) . K. Which of those constructions do you have in mind?

Download mw.mw

Worksheet: cone-map.mw

Download mm.mw

In this worksheet I call Maple's Explore() to visualize conic sections.  The graphics produced by Explore() cannot be shown on this web page.  Download the worksheet and load into Maple to play with it.

Download mw.mw

Addendum: The graphics may be improved significantly by adding a view option next to the orientation, as in:
    orientation=[115, 80, 0], view=[-3..6, -3..3, -4..4]

Maple's add() and sum() have different evaluation rules. In general, you should use add(), not sum(). if the summation range is known ahead of the time.  And if you do that, you won't need the quotation marks around a__||k.  Here is the fixed code:
restart:
p := x -> add(a__||k * x^k, k=0..5):
p(x);
p(1);
p(0);

I can't tell the purpose of the {y=-20}.  What is it supposed to do?  I am going to drop that and then fix the rest of the code.  We get:

restart;
with(PolyhedralSets):
p := PolyhedralSet([-x >= 0, -y >= 0, -z >= 0, -(1/2)*x >= 0, (-x-y)*(1/2) >= 0, (-x-z)*(1/2) >= 0]);
Plot(p);  # note the uppercase P

Download mw.mw

Download mw.mw