Does it?  The help page for ImportGraph says the supported formats are `dimacs`, `dimacs_relaxed`, `combinatorica` and `edges`.  How do these relate to graph6 and sparse6?

OK, now I think I see the reason.  Here are your equations (without the distraction of psi):

> des:= [diff(nu(xi),xi)=(theta(xi)^n - 3*nu(xi))/xi,   
   diff(theta(xi),xi)=-(nu(xi)+b*theta(xi)^(n+1))*xi*(1+b*theta(xi))/f];

Write theta(xi) = 1 + xi * u(xi), nu(xi) = 1/3 + xi * v(xi), and linearize.

> subs(theta(xi) = 1 + xi * U, nu(xi) = 1/3 + xi * V, map(rhs, des));
   R:= subs(U=u(xi), V=v(xi), map(mtaylor, %, [U=0, V=0], 2));
   L:= eval(map(lhs,des), {theta(xi) = 1 + xi * u(xi), nu(xi) = 1/3 + xi * v(xi)});
   uvdes:= zip(`=`,L,R);

uvdes := [v(xi)+xi*diff(v(xi),xi) = n*u(xi)-3*v(xi), u(xi)+xi*diff(u(xi),xi) = -(1/3+b)*xi*(1+b)/f-((1/3+b)*xi^2*b+b*(n+1)*xi^2*(1+b))/f*u(xi)-xi^2*v(xi)*(1+b)/f]

Write an equivalent second-order DE for one of the dependent variables:

> DEtools[casesplit](uvdes);

[u(xi) = -(-4*v(xi)-xi*diff(v(xi),xi))/n, diff(v(xi),`$`(xi,2)) = 1/3*(-3*n*xi*b^2-3*n*xi^2*v(xi)-n*xi-4*n*xi*b-15*n*xi^2*v(xi)*b-18*xi*f*diff(v(xi),xi)-24*v(xi)*xi^2*b^2-12*v(xi)*xi^2*b^2*n-12*v(xi)*f-16*xi^2*v(xi)*b-6*xi^3*diff(v(xi),xi)*b^2-3*xi^3*diff(v(xi),xi)*b^2*n-3*xi^3*diff(v(xi),xi)*b*n-4*xi^3*diff(v(xi),xi)*b)/xi^2/f] &where []

> vde:= op([1,2], %);

vde := diff(v(xi),`$`(xi,2)) = 1/3*(-3*n*xi*b^2-3*n*xi^2*v(xi)-n*xi-4*n*xi*b-15*n*xi^2*v(xi)*b-18*xi*f*diff(v(xi),xi)-24*v(xi)*xi^2*b^2-12*v(xi)*xi^2*b^2*n-12*v(xi)*f-16*xi^2*v(xi)*b-6*xi^3*diff(v(xi),xi)*b^2-3*xi^3*diff(v(xi),xi)*b^2*n-3*xi^3*diff(v(xi),xi)*b*n-4*xi^3*diff(v(xi),xi)*b)/xi^2/f

> dsolve(vde, v(xi), series);

v(xi) = _C1/xi*(1+(-1/30*(12*n*b+3*n+12*b+18*b^2+9*n*b^2)/f)*xi^2+1/2520*(12*n*b+3*n+12*b+18*b^2+9*n*b^2)*(18*n*b+3*n+20*b+30*b^2+15*n*b^2)/f^2*xi^4+O(xi^6))+_C2/xi^4*(12+2*(3*n*b+3*n)/f*xi^2+(-1/6*(3*n*b+3*n)*(9*n*b+3*n+8*b+12*b^2+6*n*b^2)/f^2)*xi^4+O(xi^6))+xi*((-1/30*n*(3*b^2+1+4*b)/f)+1/840*n*(3*b^2+1+4*b)/f*(6/f*n*b+10/f*b^2+5/f*b^2*n+1/f*n+20/3/f*b)*xi^2+(-1/45360*n*(3*b^2+1+4*b)/f*(6/f*n*b+10/f*b^2+5/f*b^2*n+1/f*n+20/3/f*b)*(8/f*n*b+14/f*b^2+7/f*b^2*n+1/f*n+28/3/f*b))*xi^4+O(xi^5))

Observe that the two fundamental solutions of the homogeneous equation (the  terms with _C1 and _C2) are both singular as xi -> 0 (on the order of xi^(-1) and xi^(-4)), but the particular solution of the non-homogeneous equation with _C1 = _C2 = 0 is analytic at xi=0.  So indeed it seems there is only one solution of the linearized system with your initial conditions, and it is analytic at xi=0.

To illustrate some of the subtleties that can arise here, consider a somewhat simpler linear problem:

> sys := {diff(v(t),t) = (v(t) - y(t))/t, diff(y(t),t) = v(t)};

> dsolve(sys);

{v(t) = t^(1/2)*(_C1*BesselJ(1,2*t^(1/2))+_C2*BesselY(1,2*t^(1/2))), y(t) = t*(_C1*BesselJ(2,2*t^(1/2))+_C2*BesselY(2,2*t^(1/2)))}

As usual with linear homogeneous systems, the general solution is a linear combination of two fundamental solutions. 

> sol1, sol2:= eval(%,{_C1=1,_C2=0}), eval(%, {_C2=1,_C1=0});

The one-parameter family of solutions with _C2 = 0 are analytic at t=0, i.e. can be expressed in terms of power series, with y(t) ~ _C1 * t^2/2 and v(t) ~ C_1 t as t -> 0.

> series(subs(sol1, y(t)), t);

1/2*t^2-1/6*t^3+1/48*t^4-1/720*t^5+O(t^6)

> series(subs(sol1, v(t)), t);

1*t-1/2*t^2+1/12*t^3-1/144*t^4+1/2880*t^5+O(t^6)

All other solutions have v(t) and y(t) going to a nonzero limit as t -> 0, but are not analytic at 0 because they have terms involving powers of t times ln(t).

> series(subs(sol2, y(t)), t);

(-1/Pi)+(-1/Pi)*t+(1/2*1/Pi*ln(t)-1/Pi*(-gamma+3/4))*t^2+(-1/6*1/Pi*ln(t)-1/Pi*(-17/36+1/3*gamma))*t^3+(-1/Pi*(43/576-1/24*gamma)+1/48*1/Pi*ln(t))*t^4+(-1/Pi*(-247/43200+1/360*gamma)-1/720*1/Pi*ln(t))*t^5+O(t^6)

> series(subs(sol2, v(t)), t);

(-1/Pi)+(1/Pi*ln(t)-1/Pi*(-2*gamma+1))*t+(-1/2*1/Pi*ln(t)-1/Pi*(-5/4+gamma))*t^2+(1/12*1/Pi*ln(t)-1/Pi*(5/18-1/6*gamma))*t^3+(-1/Pi*(-47/1728+1/72*gamma)-1/144*1/Pi*ln(t))*t^4+(1/2880*1/Pi*ln(t)-1/Pi*(131/86400-1/1440*gamma))*t^5+O(t^6)

You could try specifying an "initial condition" y(0) = 1, v(0) = 1, but this would not uniquely determine a solution (it would make _C2 = -Pi, but would leave _C1 undetermined).  And the solutions would all be non-analytic at 0, because of the logarithmic terms.

On the other hand, in your nonlinear problem it appears that there is an analytic solution with your initial condition, and this is unique except for the arbitrary constant psi[0].  I'm not sure why.  I'm also not sure whether there are non-analytic solutions as well, but I suspect there may be.

For example:

> eqs := [seq(eq||i, i=0..7)]; vars:= [seq(t||i, i=0..7)];

> case1:= subs(da=1,db=1,fa0=1,fa1=1,fb0=1,fb1=1,ua0=1,ua1=1,ub0=1,ub1=1,va0=1, eqs);
    solve(case1, vars);

[[t0 = -1, t1 = -1, t2 = -(-vb1*t7+vb1*vb0-vb1+vb0*t7)/(vb0*t7+vb0-1-vb1*t7), t3 = -(vb0*t7+vb0-vb1*t7-vb1)*t7/(vb0*t7+vb0-1-vb1*t7), t4 = -(vb0*t7+vb0-vb1*t7-vb1)*t7/(vb0*t7+vb0-1-vb1*t7), t5 = -(vb0*t7+vb0-vb1*t7-vb1)*t7/(vb0*t7+vb0-1-vb1*t7), t6 = (-vb1*t7^2+vb0*t7^2-vb0+1)/(vb0*t7+vb0-1-vb1*t7), t7 = t7], [t0 = -1, t1 = -1, t2 = -(2*vb0*RootOf((-2*vb1+2*vb0)*_Z^2+(-vb0+1)*_Z-vb0+1)+vb0-RootOf((-2*vb1+2*vb0)*_Z^2+(-vb0+1)*_Z-vb0+1)-1)/RootOf((-2*vb1+2*vb0)*_Z^2+(-vb0+1)*_Z-vb0+1), t3 = -RootOf((-2*vb1+2*vb0)*_Z^2+(-vb0+1)*_Z-vb0+1)-1, t4 = -RootOf((-2*vb1+2*vb0)*_Z^2+(-vb0+1)*_Z-vb0+1)-1, t5 = -RootOf((-2*vb1+2*vb0)*_Z^2+(-vb0+1)*_Z-vb0+1)-1, t6 = RootOf((-2*vb1+2*vb0)*_Z^2+(-vb0+1)*_Z-vb0+1), t7 = -RootOf((-2*vb1+2*vb0)*_Z^2+(-vb0+1)*_Z-vb0+1)]]

For example:

> eqs := [seq(eq||i, i=0..7)]; vars:= [seq(t||i, i=0..7)];

> case1:= subs(da=1,db=1,fa0=1,fa1=1,fb0=1,fb1=1,ua0=1,ua1=1,ub0=1,ub1=1,va0=1, eqs);
    solve(case1, vars);

[[t0 = -1, t1 = -1, t2 = -(-vb1*t7+vb1*vb0-vb1+vb0*t7)/(vb0*t7+vb0-1-vb1*t7), t3 = -(vb0*t7+vb0-vb1*t7-vb1)*t7/(vb0*t7+vb0-1-vb1*t7), t4 = -(vb0*t7+vb0-vb1*t7-vb1)*t7/(vb0*t7+vb0-1-vb1*t7), t5 = -(vb0*t7+vb0-vb1*t7-vb1)*t7/(vb0*t7+vb0-1-vb1*t7), t6 = (-vb1*t7^2+vb0*t7^2-vb0+1)/(vb0*t7+vb0-1-vb1*t7), t7 = t7], [t0 = -1, t1 = -1, t2 = -(2*vb0*RootOf((-2*vb1+2*vb0)*_Z^2+(-vb0+1)*_Z-vb0+1)+vb0-RootOf((-2*vb1+2*vb0)*_Z^2+(-vb0+1)*_Z-vb0+1)-1)/RootOf((-2*vb1+2*vb0)*_Z^2+(-vb0+1)*_Z-vb0+1), t3 = -RootOf((-2*vb1+2*vb0)*_Z^2+(-vb0+1)*_Z-vb0+1)-1, t4 = -RootOf((-2*vb1+2*vb0)*_Z^2+(-vb0+1)*_Z-vb0+1)-1, t5 = -RootOf((-2*vb1+2*vb0)*_Z^2+(-vb0+1)*_Z-vb0+1)-1, t6 = RootOf((-2*vb1+2*vb0)*_Z^2+(-vb0+1)*_Z-vb0+1), t7 = -RootOf((-2*vb1+2*vb0)*_Z^2+(-vb0+1)*_Z-vb0+1)]]

You first want to change to a new variable.

> f1:= simplify(eval(f, z = u*R)) assuming R>0;

     plot(f1, u=-1..1, labels=[z/R, 'f']);

First of all, your initial conditions are stated wrong.  It should be

 theta(0)=1, nu(0)=1/3

not

 theta(xi=0)=1, nu(xi=0)=1/3

But the basic problem is that dividing the right side by xi makes the system singular at xi=0.  It was not at all obvious to me that series solutions would even exist.  However, you can do something like this:

> des:= {diff(nu(xi),xi)=(theta(xi)^n - 3*nu(xi))/xi, diff(psi(xi),xi)=(n+1)*b*xi*theta(xi)^n*(1+b*theta(xi))/f,    
       diff(theta(xi),xi)=-(nu(xi)+b*theta(xi)^(n+1))*xi*(1+b*theta(xi))/f};
     N := 6;
     seriesform:= {nu(xi)=1/3 + add(nu[i]*xi^i,i=1..N),theta(xi)=1+add(theta[i]*xi^i,i=1...N),
        psi(xi)=add(psi[i]*xi^i,i=0..N)};
     map(lhs-rhs, eval(des, seriesform));
     map(series, %, xi, N+1);
     eqs:= map(t -> seq(coeff(t,xi,i),i=0..N-1),%);
     solve(eqs,indets(eqs,indexed));
     sol:= subs(%,seriesform);

sol := {nu(xi) = 1/3-1/30*n*(1+4*b+3*b^2)/f*xi^2+1/2520*n*(-20*b+60*b^3+45*b^4+70*n*b+206*n*b^2+234*n*b^3+90*n*b^4+8*n-5)/f^2*xi^4-1/408240*n*(770*b^2+840*b^3+2310*b^4+32766*b^4*n^2+24224*b^3*n^2-1566*n*b-3903*n*b^2+3780*b^5+1890*b^6+264*n*b^3+13503*n*b^4+1704*b*n^2+9186*b^2*n^2+70+17334*n*b^5+21816*b^5*n^2+5670*b^6*n^2+6615*n*b^6-183*n+420*b+122*n^2)/f^3*xi^6, psi(xi) = psi[0]+1/2*b*(n+n*b+1+b)/f*xi^2-1/24*b*(6*n*b+11*n*b^2+6*n*b^3+n^2+5*b*n^2+7*b^2*n^2+3*b^3*n^2+b+4*b^2+3*b^3+n)/f^2*xi^4+1/2160*b*(20*b^2+110*b^3+816*b^3*n^2+366*b^2*n^2-5*n-12*n*b+180*b^4+90*b^5+8*n^3+110*n*b^2+486*n*b^3+3*n^2+66*b*n^2+78*b*n^3+639*n*b^4+783*b^4*n^2+270*b^5*n^2+270*n*b^5+276*b^2*n^3+440*b^3*n^3+324*b^4*n^3+90*b^5*n^3)/f^3*xi^6, theta(xi) = 1-1/6*(1+4*b+3*b^2)/f*xi^2+1/360*(20*b+110*b^2+180*b^3+90*b^4+30*n*b+96*n*b^2+114*n*b^3+45*n*b^4+3*n)/f^2*xi^4-1/45360*(770*b^2+5600*b^3+14280*b^4+10212*b^4*n^2+7116*b^3*n^2+156*n*b+2817*n*b^2+15120*b^5+5670*b^6+12612*n*b^3+24087*n*b^4+402*b*n^2+2466*b^2*n^2+20736*n*b^5+7074*b^5*n^2+1890*b^6*n^2+6615*n*b^6-15*n+24*n^2)/f^3*xi^6}

Even in Standard, use "Maple Notation" input rather than "2D Math Notation" input.  Press Ctrl-M to switch to Maple notation for one execution group, or go to Tools, Options..., Display, switch Input display to Maple Notation, and press Apply to Session (to change it for the remainder of this Maple session) or Apply Globally to change it permanently. I rcommend the latter.

Even in Standard, use "Maple Notation" input rather than "2D Math Notation" input.  Press Ctrl-M to switch to Maple notation for one execution group, or go to Tools, Options..., Display, switch Input display to Maple Notation, and press Apply to Session (to change it for the remainder of this Maple session) or Apply Globally to change it permanently. I rcommend the latter.

Maple doesn't do the integration int(1/(x(t)+w)*y(t), w = y(t)-x(t) .. 2) because it doesn't know whether w = -x(t) is in the interval.  That would be the case if y <= 0 and  x >= -2, which is not ruled out by your assumption x(t) > y(t).  If you make the additional assumption y(t) > 0, it should work, but doesn't: I think that's because making assumptions on function calls x(t) and y(t) doesn't work as well as assumptions on variables would.  However, what you can do is use the continuous option in the int command.   Thus:

 
>  ydot := diff(y(t),t) = y(t) - int(1/(x(t)+w)*y(t), w = y(t)-x(t) .. 2, continuous);

ydot := diff(y(t),t) = y(t)+ln(y(t))*y(t)-ln(x(t)+2)*y(t)

It is unlikely that a closed form for this indefinite integral exists.  This is not necessarily a tragedy.  What would you have done next if we had given you an answer involving several pages of elliptic functions?  Practically everything you could do with that, and probably more, can be done without a closed form for the integral.  Perhaps you could use a Taylor series expansion of the integral, which should converge nicely for small x.

> J:= Int(-(-tan(0.03961111111 *Pi)+((tan(0.07800000000* Pi))^2 *x)/(sqrt(0.36-(tan(0.07800000000* Pi))^2 *x^2)))/(sqrt(1+(-tan(0.03961111111 *Pi)+((tan(0.07800000000 *Pi))^2 *x)/(sqrt(0.36-(tan(0.07800000000 *Pi))^2* x^2)))^2)),x);

> series(J, x, 20);

.1241210436*x-0.5091274650e-1*x^2-0.6535386797e-3*x^3-0.2085339117e-2*x^4-0.6306092830e-4*x^5-0.1707815782e-3*x^6-7.239812647*10^(-6)*x^7-0.1747833852e-4*x^8-9.045328537*10^(-7)*x^9-2.002915433*10^(-6)*x^10-1.188093474*10^(-7)*x^11-2.458522010*10^(-7)*x^12-1.612868265*10^(-8)*x^13-3.160641144*10^(-8)*x^14-2.241088537*10^(-9)*x^15-4.200699062*10^(-9)*x^16-3.168136559*10^(-10)*x^17-5.724666622*10^(-10)*x^18-4.538167134*10^(-11)*x^19-7.955474855*10^(-11)*x^20+O(x^21)

> plot(convert(%, polynom), x = -1 .. 1);

It is unlikely that a closed form for this indefinite integral exists.  This is not necessarily a tragedy.  What would you have done next if we had given you an answer involving several pages of elliptic functions?  Practically everything you could do with that, and probably more, can be done without a closed form for the integral.  Perhaps you could use a Taylor series expansion of the integral, which should converge nicely for small x.

> J:= Int(-(-tan(0.03961111111 *Pi)+((tan(0.07800000000* Pi))^2 *x)/(sqrt(0.36-(tan(0.07800000000* Pi))^2 *x^2)))/(sqrt(1+(-tan(0.03961111111 *Pi)+((tan(0.07800000000 *Pi))^2 *x)/(sqrt(0.36-(tan(0.07800000000 *Pi))^2* x^2)))^2)),x);

> series(J, x, 20);

.1241210436*x-0.5091274650e-1*x^2-0.6535386797e-3*x^3-0.2085339117e-2*x^4-0.6306092830e-4*x^5-0.1707815782e-3*x^6-7.239812647*10^(-6)*x^7-0.1747833852e-4*x^8-9.045328537*10^(-7)*x^9-2.002915433*10^(-6)*x^10-1.188093474*10^(-7)*x^11-2.458522010*10^(-7)*x^12-1.612868265*10^(-8)*x^13-3.160641144*10^(-8)*x^14-2.241088537*10^(-9)*x^15-4.200699062*10^(-9)*x^16-3.168136559*10^(-10)*x^17-5.724666622*10^(-10)*x^18-4.538167134*10^(-11)*x^19-7.955474855*10^(-11)*x^20+O(x^21)

> plot(convert(%, polynom), x = -1 .. 1);

This is the sort of notation typically used in thermodynamics.  You can use implicitdiff.

> implicitdiff( {q = a*x + b*y, r = b*x-a*y}, {x,y}, x, q);

a/(a^2+b^2)
  

This is the sort of notation typically used in thermodynamics.  You can use implicitdiff.

> implicitdiff( {q = a*x + b*y, r = b*x-a*y}, {x,y}, x, q);

a/(a^2+b^2)
  

There's no point in making it a function of n and m, if n and m are the dummy variables for the seq.  So just

bcs:=seq(seq(F[n,m](1) = `if`(n=0 and m=0, 1,  0), n=0..nmax), m=0..nmax),
          seq(seq(F[n,m](-1) = `if`(n=0 and m=0, -1,  0), n=0..nmax), m=0..nmax),
          seq(seq(D(F[n,m])(1) = 0, n=0..nmax), m=0..nmax),
          seq(seq(D(F[n,m])(-1) = 0, n=0..nmax), m=0..nmax);