@samiyare My comment about your replacement for infinity, which you called binfinitive, was just meant to point out that in the context its role is just like any other parameter and may cause convergence problems as the other parameters might. I didn't mean to suggest that changing variables as I did would help anything at all as far as convergence is concerned.

However, it does open up the possibility of using continuation in inf from value like 6 to N (e.g. 15) with continuation=6*(1-c)+c*N, where we replace inf with 6*(1-c)+c*N.

Something like this:

restart;
MUR:=(1-phi)^(5/2):
RhoUR:=(1-phi+phi*rho[p]/rho[f]):
RhoCPR:=(1-phi+phi*rhocp[p]/rhocp[f]):
BetaUR:=(phi*rho[p]*beta[p]+(1-phi)*rho[f]*beta[f])/(RhoUR*rho[f])/beta[f]:
dqu3:=diff(h(x),x$1)-RhoUR*BetaUR*T(x);
dqu2:=5*diff(T(x),x$2)+k[f]/k[nf]*Pr*RhoCPR*f(x)*diff(T(x),x$1);
dqu1:=5/(MUR)*diff(f(x),x$3)
+ 2*(diff(h(x),x$1)*x-h(x))
+RhoUR*(3*f(x)*diff(f(x),x$2)-diff(f(x),x$1)^2);
PDEtools:-dchange({f(x)=fy(y),T(x)=Ty(y),h(x)=hy(y),x=y*inf},[dqu1=0,dqu2=0,dqu3=0],[fy,Ty,hy,y]);
sysy:=solve(%,{diff(fy(y),y,y,y),diff(Ty(y),y,y),diff(hy(y),y)});
bcsy:={Ty(0)=1,Ty(1)=0,fy(0)=0,D(fy)(0)=inf*lambda,D(fy)(1)=0,hy(1)=0};
rho[f]:=998.2: cp[f]:=4182: k[f]:=0.597:   beta[f]:= 2.066/10000:
rho[p]:=3380: cp[p]:=773: k[p]:=36:   beta[p]:= 8.4/1000000:
k[nf]:=((k[p]+2*k[f])-2*phi*(k[f]-k[p]))/((k[p]+2*k[f])+phi*(k[f]-k[p])):
rhocp[nf]:=rho[p]*cp[p]*phi+rho[f]*cp[f]*(1-phi):
rhocp[p]:=rho[p]*cp[p]:
rhocp[f]:=rho[f]*cp[f]:
N:=15: #So here inf is really going to be 15 by continuation from 6 to 15.
sysN:=subs(inf=6*(1-c)+N*c,sysy);
bcsN:=subs(inf=6*(1-c)+N*c,bcsy);
phi:=0.0:
Pr:=7: lambda:=0.1:
res:=dsolve(sysN union bcsN,numeric,continuation=c,output=listprocedure);
[seq([y,1/N^i*diff(fy(y),[y$i])],i=0..2)];
plots:-odeplot(res,[seq([y,1/N^i*diff(fy(y),[y$i])],i=0..2)],0..1,tickmarks=[[seq(i=N*i,i=0..1,1/N)],default],labels=[x,"f, f',f''"]);
plots:-odeplot(res,[y,fy(y)-y],0..1,tickmarks=[[seq(i=N*i,i=0..1,1/N)],default],labels=[x,"f(x)-x"]);
#Alternatively:
F,F1,F2:=op(subs(res,[seq(diff(fy(y),[y$i]),i=0..2)]));
plot([F(x/N),F1(x/N)/N,F2(x/N)/10^2],x=0..15,labels=[x,"f, f', f''"]);
plot(F(x/N)-x/N,x=0..N);

@Carl Love But would you consider the non-working of Q1 a bug?

(I just submitted an SCR).
I noticed that if in Q1 the output procedure is equipped with $ then the elementwise version works:

Q11:=proc(Var::list,$) local n;
  n:=nops(Var);
  print(Var,n);
  proc(var::list,$)
       subs(Var=~[seq(1..n)],var); #Elementwise
       #subs(zip(`=`,Var,[seq(1..n)]),var);
  end proc;
end proc:
Q11([x,y]); # works

@Carl Love Thanks for pointing that out. Obviously, it should have been D(fy)(1)=0.
I have corrected it above.

In order to help we need details: the odes and the initial conditions as text or as an uploaded worksheet.

@J4James The example you give is an example for which no point computed is real. Therefore that warning and an empty plot. If on some interval lambda1 is real then plot does give a graph on that interval and nothing elsewhere: see my example above.
But if for some reason you prefer not to get the warning you could do
cdn:=alpha>=-(1/2)*S^2/(K+2);
eval(cdn,{alpha=-1,S=-0.5}); #Never satisfied for K=0..2
plot(eval(piecewise(cdn,lambda1,undefined),{alpha=-1,S=-0.5}),K=0.0..2);
but I fail to see the reason why.

@Markiyan Hirnyk I have seen you use that phrase "empty words" quite a few times. I'm not sure whether I should be offended or not. Maybe some time at your leisure you could explain what you mean by that.

But to prove the converse:
if lambda1 is a root of order >= k+1 in the polynomial p(lambda) then
p(lambda) = q(lambda)*(lambda-lambda1)^(k+1) where q(lamda) is a polynomial.
Thus (D@@i)(p)(lambda1) = 0 for i=0..k. Thus from "that equation" i.e.
 "p(D)(t^k*exp(lambda*t) ) = kth derivative of p(lambda)*(exp(lambda*t) wrt. lambda."
follows that

p(D)(t^k*exp(lambda1*t) ) = 0 for all t.

You say "My last question is: "  and yet you don't formulate a question, but repeat the assignment you have been given. I assume that you are not expecting us to do your assignment.

@Carl Love You clearly meant to write LinearAlgebra:-Eigenvectors.

You give us an image and that image doesn't even reveal how RK3 is defined.
To get reasonable help you should present the code as text or as an uploaded worksheet.

I do notice that you insist on taylor(abs(Exact-E),h=0); in order to gain insigt into the order of the error.
But recall that  g(h) = O(h^4) (as h ->0 means that there exist delta>0 and C > 0 such that

abs(g(h)) <= C*h^4  for all h satisfying abs(h) < delta.

Thus it is enough to do

taylor( Exact -E,h=0);

@Markiyan Hirnyk Here is a version which doesn't use polar coordinates. Since x' and y' both will be changed when x^2+y^2=1 and since the changes involve x' and yi for both, we need a temporary variable (temp):

res:=dsolve(sys,numeric,events=[[x(t)^2+y(t)^2-1,[temp=diff(x(t),t),diff(x(t),t)=-2*(diff(x(t),t)*x(t)+diff(y(t),t)*y(t))*x(t)+diff(x(t),t),diff(y(t),t)=-2*(temp*x(t)+diff(y(t),t)*y(t))*y(t)+diff(y(t),t)]]],range=0..200);
p1:=plot(1,theta=0..2*Pi,coords=polar,color=blue):
plots:-odeplot(res,[x(t),y(t)],0..200,axes=boxed,axes=none,refine=1):
plots:-display(%,p1);
plots:-odeplot(res,[x(t),y(t)],0..200,axes=boxed,axes=none,frames=500,numpoints=1000):
plots:-display(%,p1);

I copied the line as

Butcher, map(x->if (x=0) then `` else x end if,Butcher);

It executed without an error and returned as it should:

@Markiyan Hirnyk Sorry, I just made an addition above.

@Markiyan Hirnyk Maybe something like this:

PDEtools:-dchange({x(t)=r(t)*cos(theta(t)),y(t)=r(t)*sin(theta(t))},{diff(x(t),t,t)=0,diff(y(t),t,t)=0},[r(t),theta(t)]);
sysP:=solve(%,{diff(r(t),t,t),diff(theta(t),t,t)});
resP:=dsolve(sysP union {r(0)=.2,theta(0)=.5,D(r)(0)=1,D(theta)(0)=4},numeric,events=[[r(t)=1,diff(r(t),t)=-diff(r(t),t)]]);
p1:=plot(1,theta=0..2*Pi,coords=polar,color=blue):
p2:=plots:-odeplot(resP,[r(t)*cos(theta(t)),r(t)*sin(theta(t))],0..40,axes=none,frames=100):
plots:-display(p1,p2);

And with the range argument and maxfun=0:

resP:=dsolve(sysP union {r(0)=.2,theta(0)=.5,D(r)(0)=1,D(theta)(0)=4},numeric,events=[[r(t)=1,diff(r(t),t)=-diff(r(t),t)]],range=0..100,maxfun=0);
plots:-odeplot(resP,[r(t)*cos(theta(t)),r(t)*sin(theta(t))],0..100,refine=1,axes=none);
p1:=plot(1,theta=0..2*Pi,coords=polar,color=blue):
p2:=plots:-odeplot(resP,[r(t)*cos(theta(t)),r(t)*sin(theta(t))],0..100,axes=none,frames=500,numpoints=1000):
plots:-display(p1,p2);

@sarra Sure. But it is evident from taylor(Exact,h=0) that the first two terms are just E. Thus
E-Exact = O(h^2).

@Swhite Using NonlinearFit on your data xlist, ylist, zlist:

f:=x=a*y^b*z^c;

X:=Vector(xlist);
Y:=Vector(ylist);
Z:=Vector(zlist);
A:=Matrix([Y,Z,X],datatype=float);
res:=Statistics:-NonlinearFit(rhs(f),A,[y,z],output=solutionmodule);
res:-Results();