3. Internal code throwing an error:

solve(series(f(y)+x, x = 0, 1), y);
Error, (in series/diff) invalid arguments

4. Should this work?

solve(y^2-series(sin(x), x = 0, 2), y);
Error, (in solve/series) division by zero series

@Mac Dude

The point is to construct the inversion of the given series. This is just the simplest example. The question is why solve/series takes such a long time.

timelimit has changed in a recent version, previously it was ignored in built-in routines; but apparently now it doesn't always restore the computation state correctly (a rather obscure example is here). So perhaps it's not directly related to series.

Playing with the expansion order might be version-dependent too, but I'm getting this:

series(RootOf(G(x, y), y), x = 0, 6);
                  1-(1/2)*x^2-(1/6)*x^3+(7/48)*x^4-(1/60)*x^5+O(x^6)

series(RootOf(G(x, y), y), x = 0, 6);
Error, (in series/RootOf) unable to compute series

@vv

That's strange. Just equating the coefficients still works, with a bit of nudging:

G := (x, y) -> ln((1+x)*y)+exp(x^2+y-1)-x-cos(x);

coeff~(series(RootOf(G(x, y), y, 1), x = 0, 6), x, [$0 .. 5]) *~ [$0 .. 5]!~;
                                     [1, 0, -1, -1, 7/2, -2]

ser := eval(series(G(x, y(x)), x = 0, 6), y(0) = 1):
solve(coeff~(ser, x, [$1 .. 5]), [seq](D[1$i](y)(0), i = 1 .. 5));
   [[(D(y))(0) = 0, ((D@@2)(y))(0) = -1, ((D@@3)(y))(0) = -1, ((D@@4)(y))(0) = 7/2, ((D@@5)(y))(0) = -2]

solve/series just hangs, I'm not quite sure it's doing the right thing here. I've split it off into a new thread.

@vv

Neat. I would have probably never guessed that solve generates the inverse series (I see ?solve,series now, but it wouldn't have occurred to me to look for it). I think I discovered that subs is the way to generate the series composition pretty much by accident as well (because series uses a custom data structure, so doing substitutions on it doesn't seem particularly natural).

A couple more possible solutions:

F := (x, y) -> ln((1+x)*y)+exp(x^2*y^2)-x-cos(x):

coeff~(series(RootOf(F(x, y), y, 1), x = 0, 6), x, [$0 .. 5]) *~ [$0 .. 5]!~; # scales well

solve(coeff~(series(F(x, y(x)), x = 0, 6), x, [$0 .. 5]), [seq](D[1$i](y)(0), i = 0 .. 5)); # doesn't scale so well

plot([RootOf(F(x, y), y, 1), series(RootOf(F(x, y), y, 1), x = 0, 6)], x = -1 .. 1,
  linestyle = [solid, dot], thickness = [1, 5], view = .5 .. 1.5);

Sigh, I wish those things were more coherently documented. And also maybe that there were a forum where the developers would offer some insight on subtle issues? Apparently range does work for parametrized DEs but rather than precomputing the solution over the whole range -- like it might do in the absence of symbolic parameters -- it avoids recomputing the solution over the same interval:

f := proc(a, x) print(x); if not x::numeric then 'f(a, x)' else -a^2*sin(a*x) end if end proc:

dsol := dsolve({y(0) = 0, D(y)(0) = 1, D[1, 1](y)(x) = f(a, x)}, y(x),
  parameters = [a], events = [[diff(y(x), x) = 0, halt]], numeric, range = 0 .. Pi, maxfun = 0);

dsol(parameters = [1]);

dsol(1.);

dsol(.5);

dsol(.5) is interpolated from the already computed values. But:

dsol(2.);
Warning, extending a solution obtained using the range argument with 'maxfun' large or disabled is
highly inefficient, and may consume a great deal of memory. If this functionality is desired, it is
suggested to call dsolve without the range argument

.....print outputs.....

Warning, cannot evaluate the solution further right of 1.5707963, event #1 triggered a halt
          [x = 1.57079635131446, y(x) = .999999891130400, diff(y(x), x) = 1.90819582357449*10^(-17)]

dsol(last);
          [x = 1.57079635131446, y(x) = .999999891130401, diff(y(x), x) = 1.65762465482228*10^(-17)]

dsol(eventclear);
Error, (in _dat[1]) no events to clear

Without the range option, eventclear is actually necessary to be able to proceed.

@Parham2016

I do not follow. Your worksheet says that nu must be equal to the initial guess. So do you expect it to be equal, or close, or not necessarily close? If you do not have any assumptions about the correct value of nu, why are you dissatisfied with the value that your code generates?

And, again, your F2 isn't uniquely defined. An indefinite integral is defined up to a constant. Are you certain that it doesn't affect the result?

@ecterrab

Specifying f(z) explicitly seems to be the easiest way, thanks. When FunctionAdvisor returns a DESol object, it uses a dummy variable (not local) like _Y for the function, so I got confused here.

Can you describe your problem in more general terms? You have F1(nu, phi, r) and F3(nu, phi, r), F1 and F3 are arbitrary functions. Then you have F2(nu, phi, r) = F1(nu, phi, r)*Int(1/(xi*F1(nu, phi, xi)^2), xi = ...) -- how do you set the limits of integration here? Then you have the determinant of

 [[F1(nu, phi, 1), F2(nu, phi, 1), 0],
  [F1(nu, phi, phi), F2(nu, phi, phi), -F3(nu, phi, phi)],
  [D[3](F1)(nu, phi, phi), D[3](F2)(nu, phi, phi), -D[3](F3)(nu, phi, phi)]]

where phi is fixed. Do you want to find the values of nu where the determinant is zero? Your k is not zero at nu=1.1425.

@vv

I'm confused by this:

restart;

deq := diff(f(z), z$2)-(2*z^2+z-a-1)/z*diff(f(z), z)-((2*a+4)*z+a+1)/(2*z)*f(z):

evalf[20](DETools:-formal_sol(subs(a = .1, deq), f(z), order = 3)[1]); # correct up to the round-off
                        1.+.50000000000000000000*z+.68452380952380952380*z^2+O(z^3)

I'd understand it if in my first example Maple just returned evalf[20](cached_result). But it's doing the computation afresh and at the same time using the cached values computed at lower precision. I'm not convinced that's correct.

Come to think of it, if evalf[20](.1/2*2) extends the result beyond the original 10 digits when doing the arithmetic, why doesn't evalf[20](.1) do the same?

Maple doesn't seem to have a function to do that, something like

setDigits := (e, n) -> evalf[n](convert(e, rational, exact))

Correction: apparently `=` does work. If p and q are zppolys (not unevaluated forms like Add(zppoly, zppoly), because `=` isn't going to evaluate those), then is(p=q) gives an error the first time and FAIL when reevaluated, but evalb(p=q) seems to do the right thing, comparing p and q by content.

Also, ConvertOut(zppoly) (without a variable) does the conversion from zppoly to the coefficient list.

@acer

Yes, thanks. The issue seems to boil down to this:

proc() local x:=0 end proc(); # in 1-D Math Input
Error, reserved word `end` unexpected

proc() local x:=0; end proc(); # in 1-D Math Input
                               0

This seems to be a rather fuzzy area. I'm not even sure that local x:=0 syntax is actually documented. And if the semicolon is really a separator and not a part of a statement, perhaps the first input should work?

Moreover,

proc() local x:=0; end proc() # in 1-D Math Input

Warning, premature end of input, use <Shift> + <Enter> to avoid this message.

Supposedly the semicolon before the comment should be redundant.

@Axel Vogt

Exactly, I think that's the easiest way to see the answer.

Download inv.mw

(Does this count as using a CAS?)

@Axel Vogt

Yes. Although I meant without using a CAS.

@Axel Vogt

This might serve as a neat illustration:

compute

int(sin(t)^(1/3), t = 0 .. Pi/2) + int(arcsin(t^3), t = 0 .. 1)

in less than one minute.

@jnjn0291

Right, the eigenvalues of M won't be ordered in any nice way. But if you sort them, you get 11 smooth functions:

M := unapply(M, [x, y]):
f := proc(x, y) option remember;
  (sort@Re@LinearAlgebra:-Eigenvalues@M)(x, y)
end proc:

plot3d([seq]((i -> (x, y) -> f(x, y)[i])(i), i = 1 .. 11), 0 .. 2*Pi, 0 .. 2*Pi);

The spectrum is the set of the eigenvalues. So I don't know how you want to represent them or how you can apply implicitplot here.