Suppress the output of the loop by ending with a colon (:) instead of a semicolon (;). Unfortunately, this can only be applied to the whole loop; surprisingly, it makes no difference what punctuation you end the solve with. The colon after the loop will suppress all of its output, so now you force the output that you want by using print.

for alpha from .4 by .1 to 5 do
     S := solve(x^4-alpha, x);
     print(S[2]) ;
end do:

The smartplot is using a finer mesh or some other refinement. If you add the option gridrefine= 2 (or higher) or crossingrefine= 2 (or higher) to your implicitplot, you will get a much better plot. While these may seem like expensive options, I point out that the plot below takes 47 millisecs for me, and the plot structure only stores 104 points.

plots:-implicitplot(
     (x^2+y^2)^2 = x^3+3/10*x*y^2, x= 0..1, y= -1..1,
     gridrefine= 3, crossingrefine= 3, 
     resolution= 2^11, adaptranges,
     scaling= constrained
);

After the command you have above, put the command

print(plot1);

I think that that's all you need to do. Your command looks correct. But, not having the worksheet, I can't be sure of that.

The command ?isolve is better than msolve for this situation. But I don't think that there is any solution. Do long division. First factor 3 from the denominator:

3*p = (k^2-87)/(k+39) = k - 39 + 2*3*239/(k+39)

That leaves only a small number of possibilities for k, and I checked them all.

One of many ways to do it is with an Array-initializing procedure:

curve:= dsolve(
     {diff(y(x),x) = 1/(x-4.98), y(0)=1},
     numeric, output= listprocedure
):
(Note that I used listprocedure.)
pts:= < map(
     P-> Array(
          0..50, 
          proc(m) try P(m/10) catch: undefined end try end proc,
          datatype= float[8]
     ),
     map(rhs, curve)
)[]>^%T:

The resulting Array is plotable as is, because plot doesn't care if a few of your points are
 undefined; it just ignores them.
plot(pts);

I don't have a formula. I am guessing based on a mental analysis of the pattern. No polynomials or computerized interpolation used. The next 55 terms are

59051, 59053, 59055, 59059, 59061, 59067, 59077, 59079, 59085,
59103, 59131, 59133, 59139, 59157, 59211, 59293, 59295, 59301,
59319, 59373, 59535, 59779, 59781, 59787, 59805, 59859, 60021,
60507, 61237, 61239, 61245, 61263, 61317, 61479, 61965, 63423,
65611, 65613, 65619, 65637, 65691, 65853, 66339, 67797, 72171,
78733, 78735, 78741, 78759, 78813, 78975, 79461, 80919, 85293,
98415

You just need to save the value of omega in the same way that you save KR. I put lines of ### around the code below that I changed. Let me know if this solves your problem.

Unrelated to your problem: There's no need to use evalm. The command serves no purpose in contemporary Maple.

In the future, please post your code without the ">" beginning each line. I have to remove those one by one to use your code.

Digits:=30;

 h0:=0.156;
 d:=0.32*h0;
 l:=2;
 h1:=h0-d;
 h2:=h0+d;
 h3:=0.6*h0;
 g:=9.8;
 d1:=1;
 Term:=8;
 Num:=100:
 n:=2:
 l1:=l/n;
 
 for N from 1 to Num do
 lambda:=2*n*Pi/l: 
 k0:=evalf(10^(-10)*Pi+2.5*(N-1)*Pi/(Num-1)):
 tau0:=evalf(k0*h0):
 omega:=evalf((g*k0*tanh(k0*h0))^(1/2)):
 ###############
 Omega[N]:= omega:
 ################
 E:=g/(omega)^2:
 k1:=abs(fsolve(omega^2=g*k*tanh(k*h1),k)):
 tau1:=evalf(k1*h1):
 k2:=abs(fsolve(omega^2=g*k*tanh(k*h2),k)):
 tau2:=evalf(k2*h2):
 k3:=abs(fsolve(omega^2=g*k*tanh(k*h3),k)):
tau3:=evalf(k3*h3):
 
 u0:=evalf(g*tanh(tau0)*(1+2*tau0/(sinh(2*tau0)))/(2*k0));
 u01:=evalf(g*tanh(tau3)*(1+2*tau3/(sinh(2*tau3)))/(2*k3));
 u1:=evalf(g*(1-(tanh(tau0))^2)*(sinh(2*tau0)-2*tau0*cosh(2*tau0))/(4*(2*tau0+sinh(2*tau0))));
 H:=evalf((1+2*tau0/sinh(2*tau0))/(-lambda*k0*d));
 delta00:=evalf((lambda*d*u1/u0+I*k0)*H);
 delta01:=evalf((lambda*d*u1/u0-I*k0)*H);
 delta11:=evalf((lambda*d*u1/u0+I*k0)*H*exp(I*k0*l));
 delta12:=evalf((lambda*d*u1/u0-I*k0)*H*exp(-I*k0*l));
 delta21:=evalf(exp(I*k0*(l+l1)));
 delta22:=evalf(u0*k0*exp(I*k0*(l+l1)));
 delta31:=evalf(exp(-I*k0*(l+l1)));
 delta32:=evalf(-u0*k0*exp(-I*k0*(l+l1)));
 delta41:=evalf(exp(I*k3*(l+l1)));
 delta42:=evalf(u01*k3*exp(I*k3*(l+l1)));
 delta51:=evalf(exp(-I*k3*(l+l1)));
 delta52:=evalf(-u01*k3*exp(-I*k3*(l+l1)));
 delta61:=evalf(exp(I*k3*(l+l1+d1)));
 delta62:=evalf(u01*k3*exp(I*k3*(l+l1+d1)));
 delta71:=evalf(exp(-I*k3*(l+l1+d1)));
 delta72:=evalf(-u01*k3*exp(-I*k3*(l+l1+d1)));
 delta81:=evalf(exp(I*k0*(l+l1+d1)));
delta82:=evalf(u0*k0*exp(I*k0*(l+l1+d1))); 
 
Hi:=(Matrix([[1,1],[delta00,delta01]]))^(-1);
H1:=(Matrix([[delta21,delta31],[delta22,delta32]]))^(-1);
H2:=Matrix([[delta41,delta51],[delta42,delta52]]);
H3:=(Matrix([[delta61,delta71],[delta62,delta72]]))^(-1);
H4:=Matrix([[delta81],[delta82]]);
 
########## 
MM:= H1.H2.H3.H4:
##########
 R:=evalf(MM[2,1]/MM[1,1]):
Kr:=abs(evalf(R)):
KR[N]:=abs(Kr);
 
 end do:
 
########### L:=seq(KR[N],N=1..Num):
 
############
for N from 1 to Num do
    if KR[N]>0.2 then 
        print(Omega[N])
        ###################
    end if
end do;
 

The integral is too difficult for Maple in its current form. Convert it to spherical coordinates.

I say this with trepidation, considering how badly I erred with the limits on your previous question: Are you sure that that set R is correctly specifed? The first condition, z <= sqrt(x^2+y^2+z^2), is satisfied by any three real numbers. The second condition is strange because there is no lower limit to z. I wonder if it was supposed to be abs(z).

The integral is easily seen to be 0 without the help of Maple.

2*y >= x^2 + y^2  -->  0 >= x^2 - 2*y + y^2 = (x-y)^2  -->  (x-y)^2 = 0  -->  x=y.

Thus D is meager, and the integral is 0.

Your code contains the expression

Variance(X)=sqrt(StandardDeviation(A)).

That should be

Variance(X) = StandardDeviation(A)^2.

Also, consider using ?KernelDensity rather than EmpiricalDistribution, which gives a discrete distribution. Then you can make a custom ?Distribution by uing the function returned by KernelDensity as the ?PDF.

The two answers differ by a constant, so they are both valid forms of the integral. The constant is ln(-1) = Pi*I. You can take the derivative of each and see that in both cases you get the integrand.

It is not a silly question. The easiest way to get parameterized colours is to use colour= HUE(x) where 1 >= x >= 0. The HUE scale is the standard colour spectrum, from red=0 to violet=1.

Here's an example. I've generated example data such as you describe, but using random data since I don't have your specific vectors.

V:= unapply(Vector(16, ()-> randpoly(x, degree= 1)), x):
V||(1..4):= 'V(k)' $ 'k'= 1..4:
So now I have 4 vectors, V1, ..., V4 of 16 numeric elements each.
plot([seq]([seq]([k, (V||k)[j]], k= 1..4), j= 1..16), colour= [seq](HUE(j/16), j= 1..16));
Note that the argument to HUE is matched to the index into the Vectors.

I don't know what you mean by "both sides". Do you mean to solve it for each variable, T₀ and t? That can be done like this (I've changed the decimals to exact fractions, but it is not necessary to do that):

> 



eq:= (T[0]+exp(-sqrt(t/1000)))/T[0] = 95/100;

















> 



solve(eq, {T[0]});

















> 



solve(eq, {t});

















> 



 


















Download solve2.mw

Working from the Wikipedia pages "Gauss-Hermite quadrature" and "Hermite polynomials", I whipped up this procedure for you. It takes a positive integer argument n and returns two lists of n elements, the first being the nodes and the second being the weights.

GaussHermite:= proc(n::posint)
local 
     x,
     Hm:= unapply(numapprox:-hornerform(orthopoly[H](n-1, x), x), x), 
     R:= [fsolve](orthopoly[H](n,x), x),
     C:= evalf(2^(n-1)*n!*sqrt(Pi)/n^2)
;
     R, map(x-> C/Hm(x)^2, R)
end proc:
     
GaussHermite(6);
 [-2.35060497367449, -1.33584907401370, -0.436077411927616, 
   0.436077411927616, 1.33584907401370, 2.35060497367449], 

   [0.00453000990550894, 0.157067320322850, 0.724629595224395, 
   0.724629595224395, 0.157067320322850, 0.00453000990550894]