@J4James Do you mean something like this:

restart:
Eq1:=diff(psi(y),y$4)-diff(psi(y),y$2)=0;
res1:=dsolve(Eq1);
bcs:=psi(h1)=F,D(psi)(h1)=-1,psi(h2)=-F,D(psi)(h2)=-1;
res2:=simplify(dsolve({Eq1,bcs},psi(y)));
match(rhs(res2)=rhs(res1),y,s);
s;

@Carl Love This is a somewhat weird system. I wonder how it originated.
You cannot isolate the highest derivatives diff(u(t),t$3) and diff(w(t),t$4).
Immediately from looking at it I thought well, w and u appear differentiated to highest order 4 and 3 respectively, so you need 3+4 = 7 conditions not 6.

This extremely simple example may illustrate the point:

ode1:=diff(x(t),t)=y(t);
ode2:=diff(x(t),t,t)=diff(y(t),t);
solve({ode1,ode2},{diff(x(t),t,t),diff(y(t),t)}); #No solution
dsolve({ode1,ode2}); #Quite a few
dsolve({ode1,ode2,x(0)=1,y(0)=1}); # None it says (silently
#However
dsolve({ode1,x(0)=1});
eval(ode2,%);

@Carl Love You are right. Originally (a long time ago) the dependent variable had to be given and in the form of (example) x(t), thereby revealing both x and t. Now it is optional, but necessary in the present case but also in a case like

ode:=diff(x(t),t)=diff(h(t),t)+x(t);
dsolve({ode,x(0)=0});
dsolve({ode,x(0)=0},x(t));

I would recommend using unapply as in

g:=unapply(L(f,f)(x),x);

This has the advantage that the integration takes place at this step, whereas using

g:= x->L(f,f)(x);

no integration takes place. It takes place when a call to g is made as in g(x).

If phi(1) is negative (for some values of the parameters) and you think that your problem demands that it be positive, then either
1) your code or model is deficient
or
2) there are multiple solutions to your boundary value problem, which certainly cannot be excluded.

@samiyare I'm assuming that this is the same comment as the one I just replied to.

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