@tira I looked at your latest worksheet Exact_for_u.mw.

I corrected a few things, and did manage to get a graph in the end.
In two places a[0], which in your case is negative, appears under a square root sign. The one place is in u2, where it appears as (exp(y*sqrt(Pr*a[0])), the other place is in u7, where it appears as sqrt(-alpha[0]*a[0]).
In u2 the combined result is real (it appears). In u7, however, the result is that u7 is purely imaginary.
Thus I tried changing the sign under the square root as you will see below. Then u7 becomes real and thus also U since all the other are real.
However, there are two obvious problems: The values of u6,u7, and u8 are huge, and so is U. Also U is negative.

restart;
alpha[0] := 1/R; a[0] := (Pr-1)/(Pr*R); Gr0 := Gr/(Pr*R);
u1 := -erfc(y*sqrt(Pr)/(2*sqrt(.4)))-a[0]*((.4+(1/2)*y^2*Pr)*erfc(y*sqrt(Pr)/(2*sqrt(.4)))-y*sqrt(.4*Pr)*exp(-y^2*Pr/(4*.4))/sqrt(Pi));
u2 := (1/2)*exp(.4*a[0])*(exp(y*sqrt(Pr*a[0]))*erfc(y*sqrt(Pr)/(2*sqrt(.4))+sqrt(.4*a[0]))+exp(-y*sqrt(Pr*a[0]))*erfc(y*sqrt(Pr)/(2*sqrt(.4))-sqrt(.4*a[0])));
u3 := y*(Int(exp(-y^2/(4*R*u)+u)/(u*sqrt(u)), u = 0 .. 10))/(2*sqrt(Pi*R));
u4 := .4*a[0]*y*(Int(exp(-y^2/(4*R*u)+u)/(u*sqrt(u)), u = 0 .. 10))/(2*sqrt(Pi*R));
u5 := y*exp(.4*a[0])*(Int(exp(-y^2/(4*R*u)+u)/(u*sqrt(u)), u = 0 .. 10))/(2*sqrt(Pi*R));
u6 := y*sqrt(alpha[0])*(Int(Int(exp(-y^2/(4*R*u)+u+alpha[0]*s)*BesselI(1, 2*sqrt(u*alpha[0]*s))/(u*sqrt(s)), u = 0 .. 10), s = 0 .. 10))/(2*sqrt(Pi*R));
#u7 := y*sqrt(alpha[0]*a[0])*(Int(Int((t-s)*exp(-y^2/(4*R*u)+u+alpha[0]*s)*BesselI(1, 2*sqrt(u*alpha[0]*s))/(u*sqrt(s)), u = 0 .. 10), s = 0 .. 10))/(2*sqrt(Pi*R));
u7 := y*sqrt(-alpha[0]*a[0])*(Int(Int((t-s)*exp(-y^2/(4*R*u)+u+alpha[0]*s)*BesselI(1, 2*sqrt(u*alpha[0]*s))/(u*sqrt(s)), u = 0 .. 10), s = 0 .. 10))/(2*sqrt(Pi*R));
u8 := y*sqrt(alpha[0])*exp(.4*a[0])*(Int(Int(exp(-y^2/(4*R*u)+u+alpha[0]*s-a[0]*s)*BesselI(1, 2*sqrt(u*alpha[0]*s))/(u*sqrt(s)), u = 0 .. 10), s = 0 .. 10))/(2*sqrt(Pi*R));
Pr := .71; R := .2; Gr := .5;
evalf(eval([u1,u2,u3,u4,u5,u6],y=1.12345));
evalf(eval([u7,u8],{y=1.12345,t=1}));
U:=Gr0*(u1+u2+u3+u4-u5+u6+u7-u8)/a[0]^2;
evalf(eval(U,{y=.1,t=.1234}));
evalf(eval(U,{y=1,t=1.234}));
evalf(eval(U,{t=0,y=1}));
plot(eval(U,t=1),y=0..10);

@tira I took a look at the worksheet Exact.mw.
There are some syntax problems:

1. Many missing multiplication signs.
2. e^x is used a couple of times instead of exp(x), which is the correct syntax.
3. List brackets are used in the definition of u(y,t). They should be ordinary parentheses ().

A misprint: At the start is written alpha[o], later alpha[0].

Suggestions:

Just use the assignments  u1:= ... etc. No need for u1(y,t).
Use capital U in U:=Gr[0]*(u1+u2+u3+u4-u5+u6+u7-u8)/a[0]^2.

Most important of all: Use 1D input instead of 2D-input. It is so much easier to catch syntax errors when using 1D-math input (also known as Maple Input): In Maple go to Tools/Options/Display/Input display and choose Maple Notation. Then click the button Apply Globally (or Apply to Session).
This change applies to new worksheets. It won't change what you have written already.

Finally, there is no explanation of what is going on in the worksheet Exact.mw.. Is the claim that U solves the pde system and the boundary and initial conditions mentioned in one of your first worksheets?

@tira This is just a comment about your plots. It applies to either one of the worksheets.
I don't see the point in using seq(..., j=4), so this would be less opaque:
p[i] := pds:-plot(u, t = 4*/10, y = 0 .. 6, legend = [i]);

Another comment. You can make an animation with Gr as animation parameter like this:

Gr:='Gr': #Clearing Gr from the value assigned in the loop
pp:=proc(gr) local pds;
  pds:=pdsolve(eval(PDE,Gr=gr), IBC, numeric, spacestep = 1/100);
  pds:-plot(u, t = 4*/10, y = 0 .. 6,title=typeset("Gr = ",gr))
end proc;
plots:-animate(pp,[gr],gr=0..2,paraminfo=false,trace=24);

#You may want to leave out 'trace = 24'.

@jonlg The command

Eigenvalues(Jacobevals[2]);

should find the eigenvalues of the jacobian for the second steady state before parameters are inserted. It fails as you will notice with:

Error, (in LinearAlgebra:-Eigenvalues) expecting either Matrices of rationals, rational functions, radical functions, algebraic numbers, or algebraic functions, or Matrices of complex(numeric) values


Thus obviously this command fails too:
ev := Eigenvalues~(Jacobevals);

That command would have computed eigenvalues for all the steady states before parameters are inserted.

Thus you should insert parameters before computing eigenvalues.

You create plots one at a time, but that certainly doesn't mean that you cannot do it rather automatically.


What I would do is to write a procedure p, which takes as input a (phi, theta) pair and outputs the quantities you want corresponding to those concrete values. When that functions, i.e. when e.g. p(0.5,5); produces the quantities you want correctly, then you just have to decide how to store the results and how to go though different inputs. Using loops or seq seems easy.

eq21,eq22:=selectremove(has,eq2,Z);

@Kitonum I think the exercise consists of changing Max without using max.

I executed your worksheet in Maple 18.02 and at the end I got the result 8.998575612 from the command poleR(3, .2);

For both files the pdsolve(pde) works without problem on my machine.

@digerdiga The error message means that the procedure `DEtools/convertsys`cannot isolate the derivative of highest order (here just 1) uniquely.

This results in 2 solutions:

solve(ode,diff(x(t),t));

You may try:

debug(`DEtools/convertsys`);
DEtools[convertsys](ode,{x(0)=0},x,t);

#You may continue with examining the two different roots of solve:
undebug(`DEtools/convertsys`);
ode1,ode2:=op(op~([solve(ode,{diff(x(t),t)})]));
res1:=dsolve({ode1, x(0) = 0}, x(t), numeric);
res2:=dsolve({ode2, x(0) = 0}, x(t), numeric);
plots:-odeplot(res1,[t,x(t)],-1.1..1);
plots:-odeplot(res2,[t,x(t)],-5..5);
#####################
#You can specify the derivative at zero:
subs(t=0,x(0)=0,convert(ode,D)): %;
#Thus either D(x)(0) must be sqrt(29) or -sqrt(29).
#Now ask explicitly for a DAE method:
res:=dsolve({ode,x(0)=0,D(x)(0)=sqrt(29)},numeric,method=rkf45_dae);
plots:-odeplot(res,[t,x(t)],-1..1);

#Clearly what happens is that ode is differentiated and then the resulting second order equation is solved, as in:
odeSec:=diff(ode,t);
resSec:=dsolve({odeSec,x(0)=0,D(x)(0)=sqrt(29)},numeric);
plots:-odeplot(resSec,[t,x(t)],-1..1);

@Konstantin@ The limits of the interior integral are certainly allowed to be dependent on the integration variable of the outer integral (i.e. x).

How did you get that ymax?
ymax:=arctan(300*cos(x)+sqrt(12.25-90000*sin(x)^2));

Since x = 0..Pi the square root surely returns an imaginary result on almost the whole interval.

From looking at the result of int(f1,y) I think that your only chance for a result is numerical integration. However, you have to rethink ymax.

An example with ymax=sin(x) (nothing like yours):

Doubleint(f1,y=ymin..sin(x),x=xmin..xmax);
evalf(%);

@DimB I just get {0}:

restart;

ode1 := diff(W1(x),[x$2])+diff(gamma1(x),x)-2*g^2*W1(x) = 0;
ode2 := diff(gamma1(x),[x$2])-b*gamma1(x)-b*(diff(W1(x),x)) = 0;
ode3 := diff(Wb1(x), [x$2])+diff(gammab1(x),x)-2*gb^2*Wb1(x) = 0;
ode4 := diff(gammab1(x),[x$2])-bb*gammab1(x)-bb*diff(Wb1(x),x) = 0;

sysODE := {ode1, ode2, ode3, ode4};
solODE:=dsolve(sysODE);
odetest(solODE,sysODE);

Can you confirm that if you execute my lines above (and only those) in a new worksheet, you get as answer {0} from odetest?

That should be very simple:

subs(Sw=sum(w[k],k=1..n) ,expr);

# or

eval(expr,Sw=sum(w[k],k=1..n));

@Markiyan Hirnyk
Quote "I am somewhat troubled by no reply of you. I believe the perfect gentlemen use to answer the obtained requests."

I hate to have to repeat what I have had to do quite a few times in the past.
Your rude behavior offends me!

I regret now that I actually answered your question before this "No feedback" comment of yours.

@Markiyan Hirnyk I pointed out (what actually showed up in Kitonum's question) that before even calling simplify there was a radical (no pun intended) difference between the output from
expr and expr1 (in Kitonum's question) and from S and S1 in my simplified version.
Where in the first case (expr and S) an important and visible change happened during evaluation, in the second case (expr1 and S1) the expressions were returned seemingly untouched. In fact though, they are not returned unevaluated. sqrt has been replaced by `^`(... , 1/2) as is seen here:

type(S1,specfunc(anything,sqrt));
type(S1,`^`(anything,1/2));
type('sqrt'(4-sqrt(6)),specfunc(anything,sqrt));

My point in other words is that the procedure sqrt does some work on its own. What it does is shown in its procedure definition. I did not try to decipher what happened, but what did happen in the first case was that nested square roots were gone.