Adaptive numerical methods for DE's must have some way to prevent step sizes from getting too small, as that might cause them to run forever.  If the step size gets too small, the solver will halt with an error message.  That seems to be what has happened here.  IIRC your system had initial conditions very close to a singularity.  Since things are changing very rapidly there, the solver would need to use very small step sizes to get accurate results.  Methods such as dverk78 are very good when the equations are "nice", but may not work well close to a singularity.  Lower-order methods might be a better choice.

If you want numerical values or plots, you can do something like this:

> g := proc(y)
       if type(y, numeric) then 
          fsolve((10088-10088*exp(-x))/(72*x+76)/((1-exp(-x))/x+10088/3/(72*x+76)-1) = y, 
                     x = 0 .. infinity)
       else 'procname'(y)
       end if
    end proc;

 I don't know how you're getting that error message: it doesn't occur when I copy your input into either 1d or 2d input.  It would ordinarily mean a syntax error, with a [ occurring where it isn't allowed.  I can't figure out where that would be in your case.   In 1d input the cursor would be put at the point where this unexpected [ was found.

But I don't think you really mean that C*3[Sigma].   Yes, C*3 is C multiplied by 3.  C*3[Sigma] in 1d input is interpreted as
C * 3[Sigma] (which might not make a lot of sense, but is syntactically legal).  In 2d input it is interpreted as C*3*[Sigma].

It is very unlikely that this differential equation has any closed form solutions.  This will be the case for most nonlinear, non-autonomous differential equations that you "pull out of the air".

For the "classical" fourth-order Runge-Kutta method, you would use dsolve with the options numeric and  method=classical[rk4]);

But you'll need to give a value to the parameters q, A, m and w, or else use the parameters option.

For example:

> q:= 2: A:= 3: m:= 4: w:= 5:
  S:= dsolve([accx, accy, IC], numeric, method=classical[rk4], stepsize= 1e-14);

                      S := proc(x_classical)  ...  end;

> S(1e-10);

[t = 1.*10^(-10), x(t) = 0.383830373516018832e-3, diff(x(t), t) = 3.83831695712338435*10^6, y(t) = -9.77231644420671770*10^(-8), diff(y(t), t) = -978.032876624406982]
  

But if you don't particularly need Runge-Kutta, and there's nothing special about the system (e.g. it's not "stiff"), it's probably best to go with the default methods.

> S:= dsolve([accx, accy, IC], numeric);

                        S := proc(x_rkf45)  ...  end;

> S(1e-10);

[t=1.*10^(-10), x(t) =0.382912865089050372e-3, diff(x(t), t) = 3829141.83673023619, y(t) = -0.101523159134715600e-5, diff(y(t), t) = -10153.1532697613238]

If I understand correctly, the "real" initial conditions that your INI1 and INI2 are approximations for are
C(400/9) = 40/9, Q(400/9) = 1/80, the problem being that the right sides of both DE's would be singular there, since at
x = 400/9, c = 40/9 and q = 1/80 we have sqrt(x) - x/20 - c = 0.   To take away this singularity,
I will consider a change of independent variable, giving you an autonomous system of three DE's: write x = x(t), and note that

(dC/dt) = (dC/dx)*(dx/dt)
(dQ/dt) = (dQ/dx)*(dx/dt)

So we try:

> newsys:= {diff(x(t),t) = (sqrt(x(t)) - x(t)/20 - C(t))^2, 
diff(C(t),t) = 1/2 * (1/2/sqrt(x(t)) - 1/20 - Q(t)) * C(t)*(sqrt(x(t)) - x(t)/20 - C(t)),
diff(Q(t),t) = (Q(t) - 1/50)*(3/100 - Q(t))*(sqrt(x(t)) - x(t)/20 - C(t)) + 1/2*(Q(t)-1/40)*(1/2/x(t) - 1/20 - Q(t))*C(t)};
 

Note that for x(t) = 400/9 and C(t) = 40/9 the first two equations are automatically true, and the third becomes

diff(Q(t), t) = -(1/14400)*(40*Q(t)-1)*(31+800*Q(t))

So there is a trajectory of the new system going through the point [x=400/9. C = 40/9, Q = 1/80] along which x and C are constant while Q is not.
By the existence and uniqueness theorem for systems of differential equations, that's the only trajectory going through that point.  I conclude that there are
actually no solutions to your system approaching [C = 40/9, Q = 1/80] as x -> 400/9, whether from above or from below.  Any solution where C -> 40/9 and Q approaches some finite value as x -> 400/9 must have Q going to one of the critical points of the equation above, namely Q = 1/40 or Q = -31/800.

This sort of thing is well-studied.  For example:

Suppose G is an open set in the complex plane, z_n a sequence of points of G with no limit point in G, w_n any sequence of complex numbers.

Then there is a function A analytic in G such that A(z_n) = w_n for all n. 

Somewhat more generally (see e,g. Rudin, "Real and Complex Analysis", theorem 15.15) you can prescribe values for finitely many derivatives of f at each z_n. 

In your case with z_n = 1/n, G would not include 0.  I'm not sure what you mean by your question "Would this not uniquely define A(z) at z=0?".  Since G does not include 0, A(0) need not be defined.  And if w_n has no limit, there is no way to define A(0) that makes A continuous at 0. 

You should avoid using subscripted names as formal parameters in function definitions.   To Maple they are not just names with subscripts, they are indexed objects that could be references to components of a table or array.   If you're using 2D Math input you could right-click on all the B₁'s and C₁'s and convert them to "atomic identifiers", but that's hardly worth the trouble.  Simplest is just to call the variables B1 and C1, or B_1 and C_1 if you prefer.

A cinquefoil knot.

> with(plots):
    display([ tubeplot([3*cos(2*t)+cos(3*t), -3*sin(2*t)+sin(3*t), 1.2*sin(5*t)], t=0..2*Pi,radius=0.65, 
        tubepoints=30,numpoints=100, glossiness=0.8),
      polygonplot3d([seq([5.5*cos(2*Pi*j/5),5.5*sin(2*Pi*j/5),-2],j=0..5)],colour=black)], 
    scaling=constrained,orientation=[180,0], style=patchnogrid, lightmodel=light3);

A snail shell:

> tubeplot([t^2*sin(t),t^2*cos(t),-2*t^2], t=0..50, radius=2*sqrt(5)*Pi*t, tubepoints=40, numpoints=400, 
    scaling=constrained, orientation=[60,70],style=patchnogrid,colour=[1-t/100,t/100,1/4],lightmodel=light4);

The Kummer Quartic looks interesting.

>KQ:=-(((-1 + 3*b^2)*(1 - sqrt(2)*x - z)*(1 + sqrt(2)*x - z)*(1 - sqrt(2)*y + z)*(1 + sqrt(2)*y + z))/(3 - b^2)) + 
(-b^2 + x^2 + y^2 + z^2)^2;
 KQR := eval(KQ, {b=1.3, x=r*cos(v)*sin(w),y=r*sin(v)*sin(w),z=r*cos(w)});
 with(plots):
 implicitplot3d(KQR, r=0..3, v=0..2*Pi, w=0..Pi, grid=[50,100,100], coords=spherical, style=patchnogrid,
   lightmodel=light3,glossiness=0.8,orientation=[90,120],scaling=constrained,view=[-3..3,-3..3,-3..3]);

And, adding the option viewpoint = circleleft, we can get an animation.

I don't see those as any better than what Maple is capable of.  For example:

> plotsetup(gif,plotoptions="height=400,width=400",plotoutput="c:/test/mobius.gif");
  plot3d([cos(t)*cos((1/2)*t)*u+cos(t), sin(t)*cos((1/2)*t)*u+sin(t), sin((1/2)*t)*u], u=-2/5 .. 2/5, t = 0 .. 2*Pi, 
    grid=[4,50], scaling=constrained, orientation=[10,40], lightmodel=light3):
  plotsetup(default);

> with(plots):
  P1:=plot3d([cos(s)*(sqrt(2) *cos(s) *cos(2 *t) + 2* sin(s)* cos(t))/d(s,t),
     cos(s)*(sqrt(2)* cos(s)* sin(2 *t) - 2* sin(s)* sin(t))/d(s,t),
     3 *cos(s)^2/d(s,t)],t=0..Pi,s=-Pi/2..Pi/2, grid=[25,40], style=patchnogrid):
  P2:= plot3d([x,y,-1],x=-5..5,y=-5..5,colour=black,style=patchnogrid):
  plotsetup(gif,plotoptions="height=400,width=400",plotoutput="c:/test/boys.gif");
  display(P1,P2,scaling=constrained,orientation=[0,0],lightmodel=light4);
  plotsetup(default);

I think you want something like this:

> MyMatrix:= Matrix(N,N, (i,j) -> evalf(`if`(i=j, Diagonal(i-1), simplify(NonDiagonal(i-1,j-1)))));

The parametric form of plot3d can be used to plot a chart (or several, if the same rectangle of parameter values is used).  You can also use display in the plots package to combine several plots.

For example:
 

> P1:= plot3d( [cos(t)+1/2*cos(t)*cos(u),sin(t)+1/2*sin(t)*cos(u),1/2*sin(u)], t=0..Pi,u=0..2*Pi):
  P2:= plot3d( [cos(t)+1/2*cos(t)*cos(u),sin(t)+1/2*sin(t)*cos(u)-1,1/2*sin(u)], t=Pi..2*Pi,u=0..2*Pi):
  P3:= plot3d([[1+1/2*cos(u),v,1/2*sin(u)],[-1+1/2*cos(u),v,1/2*sin(u)]],u=0..2*Pi,v=-1..0):
  plots[display](P1,P2,P3,scaling=constrained);

> plot3d([r, theta, exp(-r)], r = 0 .. 100, theta = 0 .. 2*Pi, coords = cylindrical, axes=box);

You might want to use 10 rather than 100, because everything interesting happens for r < 10: if r is more than 6 or so you won't be able to tell the difference between z and 0 on the graph. 

Let the cards be numbered 1 to 52, so that the aces are numbers 1 to 4.  You're looking for the specific result {1,2,3,4}.

>  with(combinat):
   ti:= time(): n:= 10^6: Total:= 0:
   for count to n do
       if randcomb(52,4) = {1,2,3,4} then Total:= Total+1 end if
   end do:
   Probability = Total/n;
   Time = time()-ti;