If you are interested in calculating the area enclosed by one loop of your data set, then do the following which is a consequence of Green's Theorem in multivariable calculus.

Let (x_i,y_i), i=1,2,…,n be the coordinates of the points that form one loop of your graph.  Then the enclosed area. A, is obtained from

2A = (x₁y₂ − x₂y₁) + (x₂y₃ − x₃y₂) + … + (x_n-1y_n − x_ny_n-1) + (x_ny₁ − x₁y_n).

Change the line

diff(Q(x, t), t $ 5) = solve(ni, diff(Q(x, t), t $ 5))

to

eq := diff(Q(x, t), t $ 5) = solve(ni, diff(Q(x, t), t $ 5));

Then where you do

r := subs(diff(Q(x,t), t $5), k);

do this instead:

r := subs(eq, k);   # I think this is what you meant to do

The error goes away.  Then you are on your own for what to do with the result.

Your system of equations involves the second derivative of psi, so you need two conditions on psi, as you have prescribed.  Your system also involves the first derivatives of x and y, so you may specify x(0) and y(0), and that will yield a unique solution of the system.  You cannot arbitrarily specify extra conditions at x(1) and y(1) since in general they will conflict with the system's unique solution.

Aside: You didn't say that, but I suspect that you want to take the parameters t, r, d as unknowns which are to be determined along with the solution of the system. If that's the case, then there is a chance that a numerical solution may be found, however three parameters are too many.  You would want only two of the parameters t, r, d as unknowns, and the third one specified with a numerical value.

restart;
with(plots):
a := (k,p) -> (100-p)*exp(-k/2):
pvals := [0, 50, 100, 150, 200]:
colors := [blue, red, green, magenta, cyan]:
display(seq(pointplot([seq([k,a(k,p)], k=0..10)]), p in pvals),
    symbol=solidcircle, symbolsize=15, color=colors, legend=pvals);

Download Partial_Differential_Equations_of_Composite_Functions-with-solution.mw

Your question makes sence if x is the entire real line, and the solution and its derivatives vanish at infinity.  We cannot have real infinities in a numerical calculation, so here I am take a largish x interval and impose the vanishing conditions at the interval's endpoints.

Download mw.mw

What you have given as the "expected result" is not correct. Here is what you should get:

Download mw.mw

plot([seq(x+k, k=1..1.5, 0.1)], x=0..1, view=0..3);

dsol := dsolve({sys, ics}, numeric);
plots:-odeplot(dsol, [[t,x(t)],[t,y(t)],[t,z(t)]],
    t=0..0.005, color=[red,green,blue]);

Download mw.mw

Here is the translation into Maple of your hand-written solution.

Download mw.mw

Your command

% sh ./Maple2021.2LinuxX64Installer.run

attempts to execute a shell script but Maple2021*.run is not a shell script; it is an executable binary. So drop the "sh" and do

% chmod +x Maple2021.2LinuxX64Installer.run
% ./Maple2021.2LinuxX64Installer.run

Here is a solution along the lines of the construction that you described in response to my question. There are more sophisticated ways of doing this. This answer intentionally uses a bare minimum of Maple in order to avoid overwhelming you.

Download mw.mw

Not an equation, but a numeric procedure.

my_x := proc(t)
        if not type(t, realcons) then return 'procname'(args); end if; 
        subs(y=t, f - 0.001 = g);
        fsolve(%, x = 0 .. 1);
end proc;
plot(my_x(y), y = 0 .. 2, view = 0 .. 1, labels = [y, x]);

The equation may be solved for y and therefore the graph may be plotted through plot() rather than implicitplot(). The following produces the plot shown in Kitonum's answer.

x^2+(y-x^(2/3))^2 = 1;
solve(%, y);
subs(x^(2/3)=surd(x^2,3), [%]);
plot(%, x=-1..1, color=red, thickness=3);