## 13613 Reputation

19 years, 229 days

## You have the code...

You have the code, so what's the problem:

```restart;

ode:=diff(y(x),x)=1-2*y(x/2)^2;
ics := y(0) = 2;

y3:= dsolve({ics, ode},numeric, delaymax = 1.7);

YY5:=plots:-odeplot(y3,x=0..10);
```

## As designed...

This has been the case forever. Unassigned names are treated as nonzero by many Maple commands.
Take the simple examples:

```restart;
solve(a*x=0,x); # 0
evalb(a=0); # false
is(a=0); # false
is(a=0) assuming a=0; # true, obviously
1/a;
1/0; # error
```

For solve there is these days the parametric option:

`solve(a*x=0,x,parametric);`

In your case solve/parametric could be used like this:

```A:=Matrix([[1, 0, 0, 0, a__1], [0, 1, 0, 0, b__1], [-216, -18, 0, 0, c__1], [0, 0, 1, 0, d__1], [0, 0, 0, 1, e__1]]);
sys:=LinearAlgebra:-GenerateEquations(A,[x1,x2,x3,x4,x5],<0,0,0,0,0>);
solve(sys,[x1,x2,x3,x4,x5],parametric);
```

The answer shows correctly that the condition for the zero solution is:
1/216*c__1 + 1/12*b__1 + a__1 <> 0

## minimize and maximize...

I would begin by plotting your expression y which is clearly even, but imaginary for 1<x <2.
You don't need any readlib these days. It has been unnecessary for many years.

```restart;

y:=sqrt((x^2*(x^2-4))/(x^2-1));
###
plot(y,x=0..10);
eval(y,x=3/2);
evalc(%);
limit(y,x=infinity);
limit(y,x=1,left);
ymax:=maximize(y,x=2..10);
ymin:=minimize(y,x=2..10);
ymin:=minimize(y,x=0..1-1e-10);
```

## Singular at [x=0,y=0]...

The expression  is singular at (x,y)=(0,0):
eval(R0,{x=0,y=0}); ## error
Try another point e.g. (x,y) = (1,0):

`mtaylor(R0, [x=1, y=0], 5);`

## phaseportrait and odeplot...

You cannot get a direction field (arrows) from odeplot, but your ode is non-autonomous, so a direction field has no meaning anyway.

```restart;
ode:=diff(y(x),x,x)+8*diff(y(x),x)+9*y(x)+36*y(x)^2=2*sin(x);
## A plot of y(x) versus x:
DEtools:-phaseportrait(ode,y(x),x=-0.2..5,[[y(0)=1,D(y)(0)=0]],numpoints=1000);
# Rewriting ode as a first order system:
ode1:=diff(y(x),x)=v(x);
ode2:=subs(ode1,ode);
## A plot of diff(y(x),x) (i.e. v(x)) versus y(x):
DEtools:-phaseportrait([ode1,ode2],[y(x),v(x)],x=-.2..5,[[y(0)=1,v(0)=0]],numpoints=1000);
#### Much simpler to use dsolve/numeric directly on ode:
res:=dsolve({ode,y(0)=1,D(y)(0)=0},numeric);

## Now plot with odeplot:
plots:-odeplot(res,[y(x),diff(y(x),x)],-.2..5,numpoints=1000,thickness=2);```

## Help page...

You may have a look at the help page ?evalf,details.
Also to get inspired you could try showstat(`evalf/constant/Catalan`).
Here is a rather silly example:

```restart;
`evalf/constant/K`:=proc()  evalf(sqrt(2)*sin(sqrt(17))) end proc:
evalf(K);
evalf[15](K);
evalf[20](K);
evalf(exp(K)*Pi);
#################
showstat(`evalf/constant/Catalan`);
```

## A graphical illustration using a list...

Since I believe that Maple's solve simply uses a formula for solving polynomials of degree less than 5 the answer that comes out from solve will just be the result of inputting the coefficients of the polynomial into the formula.
Using the very same polynomial as mmcdara, here is a graphical illustration, where the returned sequence from solve is kept in a list, not a set. Thereby we are not relying on any special knowledge of the ordering of sets in Maple.
In the grapical illustration you see that none of the 3 solutions remain real: The first in the list is red, the second blue, and the third is green.

```restart;
pol:=x^3+a*x+1;
sol:=[solve](pol,x); # result: a list
with(plots):
complexplot(sol,a=-20..20,color=[red,blue,"DarkGreen"],thickness=3,style=point);
animate(complexplot,[sol,a=-20..A,color=[red,blue,"DarkGreen"],thickness=3,style=point], A=-20..20);
```

The animation follows   :

## Different...

You seem to be assuming that (a^b)^c = a^(b*c).
But that isn't correct always: Take a = -1, b=2, and c=1/3:

```((-1)^2)^(1/3);
(-1)^(2/3);
evalc(%);
```

a^b ( a <> 0) in complex analysis can be defined as exp(b*log(a)), where log is one of the branches of the logarithm. In Maple log is always chosen as the principal branch (and usually denoted by ln):

`exp(2/3*ln(-1));`

## Elementwise ~...

This isn't a procedure but it is simple:

```d:=copy(c); # if you don't want to change c
d[..,1]:=floor~(d[..,1]);
d;
c;
```

If you want to change c:

```c[..,1]:=floor~(c[..,1]);
c;```

## Using odeplot for the animation...

You get the animation using odeplot by adding frames=80:

```p1:=odeplot(Solusix,[x(t),y(t)],0..20,numpoints=1000,color=green,frames=80):
p2:=odeplot(Solusi2,[x(t),y(t)],0..20,numpoints=1000,color=red,frames=80):
p3:=odeplot(Solusi3,[x(t),y(t)],0..20,numpoints=1000,color=blue,frames=80):
display(p1,p2,p3);
```

## method=classical[rk4]...

You have an error in your ode for n(t):
You have:

```diff(n(t), t) = alpha_n*(1 - n(t)) - beta_n(t);
```

It should be:

```diff(n(t), t) = alpha_n*(1 - n(t)) - beta_n;
```

After correcting for that just do:

`dsolve(dsys1,numeric,method=classical[rk4]);`

You can specify the stepsize as in:

```res:=dsolve(dsys1,numeric,method=classical[rk4],stepsize=0.01);
res(2); #The result at t = 2```

## Syntax...

Correcting some syntactical problems I got this ODE:

`ODE:= (2/3)*int(diff(y(x),x)*x^2/(x^2 -1),x) =int(-sqrt(2*y(x)),x);`

Is this the one you intended?
If so then the answer in Maple 2023.1 (and also in Maple 2021.2 and Maple 12)  using dsolve is

`sqrt(y(x)) + 3/4*sqrt(2)/x + 3/4*sqrt(2)*x - c__1 = 0`

Differentiating ODE you get a separable ode:

`ode2:=diff(ODE,x);`

From which the resulting implicit solution follows by separation of variables and then integrating both sides wrt. x.

But my interpretation of your original ODE may be wrong?

## indets/float...

To get the floats do:

`indets(X1,float);`

## plottools:-getdata...

Save the plot itself in some variable. Below I use p.

```p:=%; # saving the plot right after executing the plot3d command.
data:=plottools:-getdata(p);
data[1];
data[2];
data[3];
dx:=5./48; #spacing between x-values
dy:=evalf(b)/48; ##spacing between y-values
X:=[seq(dx*i,i=0..48)]; # List of x-values
Y:=[seq(dy*i,i=0..48)]; # List of y-values
[X[7],Y[9]]; # An example xy point:
data[3][7,9]; # The corresponding value of phi[2]
eval(phi[2],{x=X[7],y=Y[9]}); # Check
```

## Telescoping...

As you notice this series is telescoping. So this works:

```restart;
A:=sum(1/f(n)-1/f(n+1),n=1..infinity); # Using the unassigned f instead of sqrt
B:=subs(infinity=N,A);
expand(B);
eval(%,f=sqrt);
limit(%,N=infinity);
```

 2 3 4 5 6 7 8 Last Page 4 of 160
﻿