One way is to use fdiff and use output=listprocedure like this.

res:=dsolve(numeric,procedure=p,initial=Array([1,0,0,1]),number=4,procvars=[x(t),diff(x(t),t),y(t),diff(y(t),t)],start=0,output=listprocedure);
#As before:
plots:-odeplot(res,[[t,x(t)],[t,y(t)]],0..5);
#Now pull out the separate procedures for x(t), y(t), diff(x(t),t), and diff(y(t),t).
X,Y,X1,Y1:=op(subs(res,[x(t),y(t),diff(x(t),t),diff(y(t),t)]));
#Then use fdiff:
?fdiff

plot(fdiff(X1,[1],t),t=0..5);

One way is to use fdiff and use output=listprocedure like this.

res:=dsolve(numeric,procedure=p,initial=Array([1,0,0,1]),number=4,procvars=[x(t),diff(x(t),t),y(t),diff(y(t),t)],start=0,output=listprocedure);
#As before:
plots:-odeplot(res,[[t,x(t)],[t,y(t)]],0..5);
#Now pull out the separate procedures for x(t), y(t), diff(x(t),t), and diff(y(t),t).
X,Y,X1,Y1:=op(subs(res,[x(t),y(t),diff(x(t),t),diff(y(t),t)]));
#Then use fdiff:
?fdiff

plot(fdiff(X1,[1],t),t=0..5);

The number of elements in R is the number of function evaluations stored in the remember table of the procedure q. This will contain the information that input 't' returns literally q(t). It will have 2 entries for zero, one for the exact 0 and one for the float 0. (unless y(0)= 1 is replaced by y(0.) = 1).
My code attempted to find the actual stepsize used, so I removed the entry for 't'.

To find the number of steps for a non-classical method from the number of function evaluations, you have to know about the workings of the particular method.

nops(L) gives the number of operands in L (=elements, when L is a list as I used). I didn't use a vector, but if you want that you could just change the definition of Lix to

Lix:=<remove(type,map(op,[indices(R)]),name)>:

or equivalently

Lix:=Vector(remove(type,map(op,[indices(R)]),name)):

In Maple 15 the number of elements in a vector (or list, set, table, matrix) can be found by

numelems(Lix);

numelems(R);

In earlier versions you could e.g. do

LinearAlgebra:-Dimension(Lix);
nops({indices(R)});

Notice that the package linalg and the corresponding structures vector, matrix, array are deprecated. You should use the package LinearAlgebra and Vector, Matrix, and Array instead.

The number of elements in R is the number of function evaluations stored in the remember table of the procedure q. This will contain the information that input 't' returns literally q(t). It will have 2 entries for zero, one for the exact 0 and one for the float 0. (unless y(0)= 1 is replaced by y(0.) = 1).
My code attempted to find the actual stepsize used, so I removed the entry for 't'.

To find the number of steps for a non-classical method from the number of function evaluations, you have to know about the workings of the particular method.

nops(L) gives the number of operands in L (=elements, when L is a list as I used). I didn't use a vector, but if you want that you could just change the definition of Lix to

Lix:=<remove(type,map(op,[indices(R)]),name)>:

or equivalently

Lix:=Vector(remove(type,map(op,[indices(R)]),name)):

In Maple 15 the number of elements in a vector (or list, set, table, matrix) can be found by

numelems(Lix);

numelems(R);

In earlier versions you could e.g. do

LinearAlgebra:-Dimension(Lix);
nops({indices(R)});

Notice that the package linalg and the corresponding structures vector, matrix, array are deprecated. You should use the package LinearAlgebra and Vector, Matrix, and Array instead.

There is still a ws inIBCZ in your corrected version.

There is still a ws inIBCZ in your corrected version.

Try raising Digits above 10, then fsolve finds nothing.

restart;
Digits:=11:
fsolve([2.32=a+b,pa=.4310344828*a,pb=.4310344828*b,pa+pb=1]);

(See perhaps the addendum to my answer 'No solution indeed').

Try raising Digits above 10, then fsolve finds nothing.

restart;
Digits:=11:
fsolve([2.32=a+b,pa=.4310344828*a,pb=.4310344828*b,pa+pb=1]);

(See perhaps the addendum to my answer 'No solution indeed').

@icegood Making rather obvious changes in the code for IntegrationTools:-Expand we get the following procedure ExpandSum. To see the code for IntegrationTools:-Expand, do showstat(IntegrationTools:-Expand);

ExpandSum := proc(v)
local f, i, c, g, w;
      if type(v,specfunc(anything,{'Sum','sum'})) then
        f := expand(op(1,v));
        i := op([2,1],v);
        if type(f,`*`) then
          g, c := selectremove(type,f,(':-dependent')(i));
          c*op(0,v)(g,op(2 .. -1,v))
         elif type(f,`+`) then
          return map(procname,map(op(0,v),f,op(2 .. -1,v)))
         else
          return v
         end if
       else
        if not hastype(v,specfunc(anything,{'Sum','sum'})) then
          return v
         else
          return evalindets(v,specfunc(anything,{'Sum','sum'}),'procname')
        end if
       end if
end proc;

@icegood Making rather obvious changes in the code for IntegrationTools:-Expand we get the following procedure ExpandSum. To see the code for IntegrationTools:-Expand, do showstat(IntegrationTools:-Expand);

ExpandSum := proc(v)
local f, i, c, g, w;
      if type(v,specfunc(anything,{'Sum','sum'})) then
        f := expand(op(1,v));
        i := op([2,1],v);
        if type(f,`*`) then
          g, c := selectremove(type,f,(':-dependent')(i));
          c*op(0,v)(g,op(2 .. -1,v))
         elif type(f,`+`) then
          return map(procname,map(op(0,v),f,op(2 .. -1,v)))
         else
          return v
         end if
       else
        if not hastype(v,specfunc(anything,{'Sum','sum'})) then
          return v
         else
          return evalindets(v,specfunc(anything,{'Sum','sum'}),'procname')
        end if
       end if
end proc;

If you want to use lambda as a parameter you somehow have to avoid it appearing as an argument to f as it does in one of the boundary points.

So replace x with x(t) satisfying x'(t) = 1 and x(0) = -lambda, and let g(t) = f(x(t)). Then you get the following system

sys:={g(t)^2*diff(g(t),t,t)+x(t)*g(t)*(1-2*beta)+beta*diff(g(t),t)*(x(t)^2-1)=0,diff(x(t),t)=1};
#Now we use 3 parameters p (as in my previous answer), lambda, and beta.
#Otherwise the idea is as before although there are some changes.
init:=dsolve(sys union {g(0)=beta*(lambda^2-1)/p,D(g)(0)=p,x(0)=-lambda},numeric,parameters=[p,lambda,beta],output=listprocedure);
G,X:=op(subs(init,[g(t),x(t)]));
#We don't really need X, but use it as a check.
q:=proc(P,L,B)
   if not type([P,L,B],list(numeric)) then return 'procname(_passed)' end if;
   init(parameters=[P,L,B]);  
   G(L+1);
end proc;
q(.0236,.3,.01);
#I use the plotting data to find the value of p where G(lambda+1)=0 (i.e. f(1)=0).
plot(q(p,.3,.01),p=.023..0.024,-1..1);
op(indets(%,listlist)):
L:=remove(has,%,HFloat(undefined)):
Ls:=sort(L,(x,y)->evalb(abs(x[2])Ls[1];
p1,q1:=op(Ls[1]);
# fsolve has problems with fsolve(q(p,.3,.01),p=p1), which is not so surprising since f'(1) is infinite when f(1) = 0. So the following is just a check.
fsolve(q(p,.3,.01)-q1,p=p1);
X(0);
G(0);
G(.3+1);
#Another value for lambda, same beta.
plot(q(p,.35,.01),p=0..0.2,-1..1);
op(indets(%,listlist)):
L:=remove(has,%,HFloat(undefined)):
Ls:=sort(L,(x,y)->evalb(abs(x[2])Ls[1];
p1,q1:=op(Ls[1]);
fsolve(q(p,.35,.01)-q1,p=p1);
X(0);
G(0);
G(.35+1);

This method could be automated, but maybe somebody has a better idea.

Added: Instead of the plotting idea which makes use of adaptive plotting you could use the following procedure R, which takes lambda, beta, and a p-interval as its first 3 arguments and outputs a p-value with abs(f(1))< a tolerance, by default 10^(-Digits+2).

R:=proc(L::numeric,B::numeric,interval::range(numeric),{initstep::positive:=1e-2,tolerance::positive:=10^(-Digits+2),iterationlimit::posint:=10})
  local i,j,pp,M,a,b,d,A;
  a,b:=op(interval);
  d:=initstep;
  A:=tolerance+1;
  for j from 1 to iterationlimit while A>tolerance do
    M:=Matrix(0);
    i:=0;
    for pp from a to b by d do
      try
        i:=i+1;
        M(1,i):=pp;
        M(2,i):=q(pp,L,B);
      catch:
        break
      end try
    end do;
    if i=1 then error "Try another interval" end if;
    A:=abs(M(-3));
    a:=M(-4);
    d:=d/10;
  end do;
  if A>tolerance then WARNING("Result not within tolerance") end if;
  M(-4),M(-3);
end proc;
R(.35,.01,0.01..1);
R(.36,.01,0.01..1);
# A loop through lambda-values, initially crude since f(1) may never become zero.
for L from .3 to .5 by .02 do L,R(L,.01,.01..1,iterationlimit=2) end do;
for L from .3 to .37 by .01 do R(L,.01,.01..1) end do;


R(.35,.2,.01..1,iterationlimit=7,tolerance=1e-3,initstep=.01);
R(.35,.5,.01..1,iterationlimit=7,tolerance=1e-3,initstep=.07);

If you want to use lambda as a parameter you somehow have to avoid it appearing as an argument to f as it does in one of the boundary points.

So replace x with x(t) satisfying x'(t) = 1 and x(0) = -lambda, and let g(t) = f(x(t)). Then you get the following system

sys:={g(t)^2*diff(g(t),t,t)+x(t)*g(t)*(1-2*beta)+beta*diff(g(t),t)*(x(t)^2-1)=0,diff(x(t),t)=1};
#Now we use 3 parameters p (as in my previous answer), lambda, and beta.
#Otherwise the idea is as before although there are some changes.
init:=dsolve(sys union {g(0)=beta*(lambda^2-1)/p,D(g)(0)=p,x(0)=-lambda},numeric,parameters=[p,lambda,beta],output=listprocedure);
G,X:=op(subs(init,[g(t),x(t)]));
#We don't really need X, but use it as a check.
q:=proc(P,L,B)
   if not type([P,L,B],list(numeric)) then return 'procname(_passed)' end if;
   init(parameters=[P,L,B]);  
   G(L+1);
end proc;
q(.0236,.3,.01);
#I use the plotting data to find the value of p where G(lambda+1)=0 (i.e. f(1)=0).
plot(q(p,.3,.01),p=.023..0.024,-1..1);
op(indets(%,listlist)):
L:=remove(has,%,HFloat(undefined)):
Ls:=sort(L,(x,y)->evalb(abs(x[2])Ls[1];
p1,q1:=op(Ls[1]);
# fsolve has problems with fsolve(q(p,.3,.01),p=p1), which is not so surprising since f'(1) is infinite when f(1) = 0. So the following is just a check.
fsolve(q(p,.3,.01)-q1,p=p1);
X(0);
G(0);
G(.3+1);
#Another value for lambda, same beta.
plot(q(p,.35,.01),p=0..0.2,-1..1);
op(indets(%,listlist)):
L:=remove(has,%,HFloat(undefined)):
Ls:=sort(L,(x,y)->evalb(abs(x[2])Ls[1];
p1,q1:=op(Ls[1]);
fsolve(q(p,.35,.01)-q1,p=p1);
X(0);
G(0);
G(.35+1);

This method could be automated, but maybe somebody has a better idea.

Added: Instead of the plotting idea which makes use of adaptive plotting you could use the following procedure R, which takes lambda, beta, and a p-interval as its first 3 arguments and outputs a p-value with abs(f(1))< a tolerance, by default 10^(-Digits+2).

R:=proc(L::numeric,B::numeric,interval::range(numeric),{initstep::positive:=1e-2,tolerance::positive:=10^(-Digits+2),iterationlimit::posint:=10})
  local i,j,pp,M,a,b,d,A;
  a,b:=op(interval);
  d:=initstep;
  A:=tolerance+1;
  for j from 1 to iterationlimit while A>tolerance do
    M:=Matrix(0);
    i:=0;
    for pp from a to b by d do
      try
        i:=i+1;
        M(1,i):=pp;
        M(2,i):=q(pp,L,B);
      catch:
        break
      end try
    end do;
    if i=1 then error "Try another interval" end if;
    A:=abs(M(-3));
    a:=M(-4);
    d:=d/10;
  end do;
  if A>tolerance then WARNING("Result not within tolerance") end if;
  M(-4),M(-3);
end proc;
R(.35,.01,0.01..1);
R(.36,.01,0.01..1);
# A loop through lambda-values, initially crude since f(1) may never become zero.
for L from .3 to .5 by .02 do L,R(L,.01,.01..1,iterationlimit=2) end do;
for L from .3 to .37 by .01 do R(L,.01,.01..1) end do;


R(.35,.2,.01..1,iterationlimit=7,tolerance=1e-3,initstep=.01);
R(.35,.5,.01..1,iterationlimit=7,tolerance=1e-3,initstep=.07);

No problem with your last attempt except of course if i doesn't have a numerical value, or if k does have a numerical value.

Example.

restart;
i:=5:
for j to i do c[j]:=plot(x^j,x=0..1) end do:
plots:-display(c[k]$k=1..i);

#Since j has a numerical value (in this case 6) you need unevaluation quotes as shown below, therefore the use of seq in place of $ (as suggested by Clare So) is in general better.
plots:-display('c[j]'$'j'=1..i);

The elementwise operation ~ used in Robert Israel's answer was introduced in Maple 13. You can use map instead.

plots[display](map(plottools[scale],{p1,p2,p3},1/(2*Pi*sqrt(2)),1),
            labels=[x/(2*Pi*'sqrt(2)'),h]);

The elementwise operation ~ used in Robert Israel's answer was introduced in Maple 13. You can use map instead.

plots[display](map(plottools[scale],{p1,p2,p3},1/(2*Pi*sqrt(2)),1),
            labels=[x/(2*Pi*'sqrt(2)'),h]);