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

Is there a reason for looking at xhat as a function of all those variables?  Why not take xhat to be a function of a_1 only.  The solution will involve the remaining variables as parameters.

Download mw.mw

There is no meaningful way to do what you are asking, and that's not because of any shortcoming on the part of Maple.

Your equations are defined in a three-dimensional phase space, that is, the associated vector field is 3D and the solution curve is a curve in 3D.  You plot the projection of that 3D curve onto the (V,In) plane.  That works fine since the projection of a 3D curve onto a 2D plane is still a curve.

In contrast, the projection of a 3D vector field onto a 2D plane is not a vector field because each point of your 2D plotting domain is the foot of the projection of all the vectors in 3D that lie "above" that point.

Change
if not is({args}, set(numeric)) then return ('procname')('args') end if;
to
if not type({args}, set(numeric)) then return ('procname')('args') end if;
and it will work.

Comment added later:

It is likely that numeric is not what you want.   For instance, Pi is not numeric.  To allow for symbolic constants such as Pi as arguments, change numeric to realcons.  Here is the way I would do it:

if not type([args], list(realcons)) then return 'procname'(args) end if;