You can pass to dsolve an Array of values of the independent variable. Then dsolve's output includes a Matrix whose columns are the independent and dependent variable values.

Sol:= dsolve(
     {diff(y(t),t$2) - y(t)^2 = 1, y(0)=0, D(y)(0)=0},
     {y(t)},
     numeric,
     output= Array(1..101, k-> .01*(k-1))
)[2,1][.., [1,2]]:
 
ExcelTools:-Export(Sol, "filename"):

See ?ExcelTools,Export for a few options on that command.

If you are sure that you want to turn it into a fraction, then use convert(..., rational) or convert(..., rational, exact). But Maple can also sometimes identify an irrational expression that yields a given decimal. Such is the case with the decimal that you gave.

identify(0.569840290998053);

Sorry that I missed the coefficient 2 before. With that, there are no rational roots. So the exact answer is substantially more complicated, but the steps to getting that answer are essentially the same. You'll need to download the attached worksheet to see the full answer, because it won't display in MaplePrimes.

Download parfrac.mw

I don't know what went wrong with the factor command that you gave. When I give it, I get

f:= x-> x^3+x^2+x+1:
factor(f(x));

I suggest that you try the factor again, after doing a restart. I also suggest that you explicitly type the command rather than cut-and-paste it. You also must have made a mistake with your synthetic division because -1 is a root.

Continuing, we note that the numerator and denominator have a common factor x+1. After this reduction, it is no longer a partial-fractions problem.

The answer is 4, and that should be obvious to you. But we can nonetheless simulate this situation with Maple.

#Generate random augmented matrix.
A:= LinearAlgebra:-RandomMatrix(10,15):
R:= LinearAlgebra:-GaussianElimination(A, method= FractionFree):
#Verify last row is not all zero or inconsistent:
R[10, 14..15];

Sol:= LinearAlgebra:-LinearSolve(R, free= t):
#Count free variables in solution.
nops(indets(Sol, name));
    
     4

I think that you are saying that you want the various orderings of the compositions to be considered as one for the purpose of random selection. Then you want to use partition instead of composition, like this:

RandomPartitions:= proc(n::posint, k::posint)
local
     C,
     Partitions:= [seq(C-~1, C= select(p-> nops(p)=k, combinat:-partition(n+k,n+1)))],
     Rand:= rand(1..nops(Partitions))
;
     ()-> Partitions[Rand()]
end proc:

R:= RandomPartitions(8,6):
n:= 1000:
map(lhs=rhs/n, Statistics:-Tally(['R()[]' $ n]));

Unfortunately, "approximate using n digits precision" means that the subcomputations are done to n digits, not that the overall computation is accurate to n digits. It is the responsibility of the user to increase the precision of the subcomputations if that is needed.

Your code is selecting the second point of the two returned by solve. Apparently, you want the first point. So, you need a criterion with which to select between the two points returned by solve. Would the point which is closest to the origin be an appropriate criterion? Then you can use this procedure.

Closest:= proc(
     Sol::list(list(name=numeric)),
     {To::list(numeric):= [0,0],
      XY::list(name):= ['x', 'y']
})
local
     Min:= infinity, ptMin, pt, d,
     dist:= proc(Pt1,Pt2)
     option inline;
          evalf(sqrt(`+`(((Pt1 -~ Pt2)^~2)[])))
     end proc
;
     for pt in Sol do
          d:= dist(eval(XY, pt), To);
          if d < Min then
               Min:= d;
               ptMin:= pt
          end if
     end do;
     
     ptMin
end proc: 

Pt:= Closest(Sol);

Point2:= plottools:-disk(eval([x,y], Pt), 0.2e-1, color = red):

How about this?

n:= 10:
map(lhs=rhs/n, Tally(['R()[]' $ n]));

seq(coeff(eq,v), v= vars);

Using ex as your original expression,

subsop(1= factor(op(1, ex)), ex);

convert(expand(u^3+u^2-1), exp);

There are numerous problems in your code. First, the inverse tangent function is arctan, not atan. Second, you cannot make an assignment inside a return statement (or any other statement). Unlike C, in Maple an assignment is not an expression that can be part of another expression. Third, what's the point of having a loop if you're always going to return on the first iteration?

Correct these and then maybe we can discuss the rest.

I was able to get a result by substantially reducing the requested accuracy. I had to increase abserr, but not relerr. So, I'm a bit skeptical about this result. I changed your theta to T to save typing.

restart:
ODE:=
     (1+T(x))*diff(T(x),x) + x*diff(T(x),x)^2 +
     x*(1+T(x))*diff(T(x),x$2) - T(x)^3-diff(T(x), x)
     = 0
;

BCs:= T(1)=1, D(T)(0) = 0:
Sol:= dsolve(
     {ODE, BCs}, numeric, method= bvp[middefer],
     maxmesh= 2^15, abserr= .58e-1
):
plots:-odeplot(Sol, [x,T(x)], x= 0..1);

subs is a low-level command that can only be used to substitute for subexpressions that appear as nodes on the tree that represents the expression. So using subs requires some (small) amount of knowledge of how expressions are stored. For example, in the expression x+y+z, the nodes are x+y+z, x, y, and z; x+y is not a node. You can see the tree representing an expression by using ToInert(expression). Unfortunately, the representation is not graphical.

algsubs is a high-level command that can be used to substitute for more-complicated subexpressions. The syntax is the same as subs. In your case, that would be

algsubs(cosh(k*h)/sinh(k*h)= 1/tanh(k*h), td);