Tried your code snipruns  in Maple 18, Maple 2015, Maple 2016, Maple 2017, and Maple 2018. No errors in any of them.

'mverbatim' is something which seems to crop up when Maple is unable to parse 2-D input.

From the above, I have only have 2 suggestions

1. Post a more complete piece of code, beginning with a restart() command which exhibits the problem
2. Don't use 2-D input :-)

Type

?Maple Programming Guide

at the Maple prompt. This should open a new worksheet.

In this worksheet, you want the "Advanced Connectivity" section (section 14, for current Maple installations). Then check out the subsection titled "The Maple Command Line Interface" (subsection 14.4, for current Maple installations)

I can only suggest a careful reading of the pdsolve/numeric help page. At the risk of oversimplifying this "help" page

1. The pdsolve(...,numeric,..) command will return a 'module', comprising four "methods" for displaying/manipulating solutions. These "methods" are
   1. plot3d(). The pdsolve(.., numeric,..) command is restricted to systems of 2 (or fewer) independent variables. With 2 independent variables, plotting one or more of the dependent variables of the system requires a 3-D plot, which is what you get with this method. You can also plot (more-or-less any) function of the dependent and independent vaiables
   2. plot(): This more or less produces 'slices' through the 3-D plots produced above: ie fix one of the independent variables, then plot some function of the dependent and independent variables, whilst varying the other independent variable
   3. value() for producing numerical values
   4. animate() for producing animations (obviously)
2. The attached merely scratches the surface of the possibilities for your specific problem.

Download pdesol.mw

as in the attached

Download soleq.mws

You mean something like?

Download morph.mw

Wee I do a simple plot of your equation over a very restricted range of values for x, -0.01..0.01, and I get the figure shwn in the attached. If the *solution* fo your equation is the point at which the dependent variable is zero - then you have an awful lot of these, even in thi ting range of of the independent variable. Before I go any further, please define a lot more carefully what you mean by the statement

" How can I solve this equation?"

I mean you want all of the values where the function is zero? some of them?, any one in particular?

Download cplot.mw

If I just take the code you posted, and your requirement that " I can't get all of the values of r out of the loop ", then I can get all the 'r-values' out of the loop. This may not be the best way to solve your problem, but I had to make so many guesses about the intent of your code, that it's the all I'm prepared to go with right now.

Note with all the guesses and asssumtions I had to make, 'r' comes out as 8, every time: change to the value of 't' and the definition of the procedure cumprob() may change this - but since you do not define these I  just have to guess.

#
# Need some random value for 't' just to make this
# worksheet executable - so I picked t-3
#
   t:=3:
#
# Need some random procedure called cumprob() just to
# make this worksheet executable - so I picked a simple
# squaring function
#
   cumprob:= x-> x^2:
#
# Load the NumberTheory package, just to make this
# worksheet executable
#
   with(NumberTheory):
#
# Now everyth8ing at least "functions" - wtf was the OP's
# problem again?? Seems to want all the values of 'r'
# which get generated in this loop: seems not to realise
# that every time (s)he computes 'r', it will overwrite
# the previous value. So specify 'r' as NULL, which just
# initialises a sequence then for each computed value of
# 'r' , add it to the *sequence* of values in 'r'
#
  r:=NULL:
  for i from 1 to 10 do
        realchoose := rand(0.0..1.0);
        o := realchoose();
        print(o);
        for q from 1 to 2^t do # this won't execute cos 't' is undefined
              if      cumprob(q)>=o # this won't execute cos cumprob is undefined
              then printf("measurement = %d\n", q-1); break;
              end if;
        end do;
        if      q-1=0
        then q:=2^t;
        end if;
        p := ContinuedFraction((q-1)/2^t);# this won't work unless the NumberTheory package is loaded
        cf := Term(p, 1..10);
        if       cf=0
        then cf:=1;
        end if;
        print(cf);
        cfmax := max[index](cf);
        cfmax1 := max(cf);
        if      cfmax1<10
        then cfmax:=5;
        end if;
        R := Convergent( p, cfmax-1);
        r := r, denom(R); # append the current value of denom(r) to previous
                                     # values
       print(r);
   end do:
#
# Now check the sequence of values obtained for 'r'
#
  r;

Maple is quite 'polite' about writing files in certain locations - after all you woulldn't want Maple to overwrite your whole "C" drive.

You may need to run Maple with administrator privileges in order to (over)write your Maple installation. However you normally start Maple (with a left-click), use a right-click and select "Run as Administrator" - then be very, very careful.

For some of the values of rho which you specify, the argument for the second boundary condition is complex.

In general you cannot solve complex-valued BVPs in Maple. You can solve complex-valued IVPs in Maple, so a shooting method construction might be possible.  However, I thought (I might be wrong) that even Maple's IVP solver was restricted to the case where the independent valiable is real, and one or more of the dependent variables is complex

One can trick the BVP solver for cases where only the dependent variable is complex - just by redefining it as a sum of real and imaginary parts, thus generating two real equations which can be solved as a system. However, I'm pretty sure that there is no equivalent approach for when the independent variable is complex.

Note that in 100000 iterations, one gets values in ranges 0..1, all the way up to 16..17. The following counts all of these

restart;
f:=x->1-2*(1/x-trunc(1/x))+1/x;
x:=50100;
A:=Array(0..20):
for j from 1 to 100000 do
      tmp:=f(x);
      p:=floor(evalf(tmp));
      A[p]:=A[p]+1;
      x:=tmp
od:
A;

makes it more understandable

Download crits.mw

if you upload the worksheet you have using the big green up-arrow in the Mapleprimes toolbar.

You also might want to check the help for the command

VectorCalculus:-Flux()

just type ?VectorCalculus:-Flux at the Maple prompt. As the help states, this will

compute the flux of a vector field through a surface in R^3 or a curve in R^2

restart;
BigProc:= proc(H::list,T::list)
                         local Form, i;
                         Form:=NULL;
                         for i from 1 to nops(H) do
                              Form := Form,(1 + add(H[j]*T[j],j=1..i))/(1 + (add(1/(H[j]*T[j]),j=1..i)));
                         end do:
                         return [Form];
                end proc:
BigProc([1,3,5],[3,6,8])

depending on the form you want.

restart;
solution := {p[1] = 2.788944999, p[2] = 4.940143518}, { p[1] = 15.29764736, p[2] = 4.946617373};
M1:=Matrix(convert~([solution], list));
M2:=rhs~(Matrix((convert~([solution], list))));

f:=piecewise(x<1, 2*x,
                       x=1, 4,
                       x>1, 3*x
                     );
plot( f, x=0..2, discont=true);