Another comment.

If you do

solve(RedSys,{diff(X74(t),t,t,t),diff(X97(t),t,t),diff(X63(t),t,t),diff(X21(t),t),diff(X16(t),t)});
 

NULL is returned.

Thus Maple cannot solve for the leading derivatives, here meaning the highest derivatives of the 5 functions appearing in the 5 equations.

Preben Alsholm

If you want to run muller.mw (not mulleralg.mw) as it is written, you need the library 'nalib',  which is available at the top of the page
http://www.cbu.edu/~wschrein/pages/M329Notes.html

Choose 'self-extracting Windows archive'. Save it in a convenient location. You need to know the path. Say it is

"F:/mymath/nalib"

then in Maple you do

restart;
libname:="F:/mymath/nalib", libname;
with(numanal);
#then you are ready.
Make sure that in the path you use `/` or  `\\`, not `\`.
 

The worksheet mulleralg.mw doesn't need the library because the procedure definition is right in there.

Walter Schreiner's numanal package contains many other procedures too.

Preben Alsholm

You may experiment with something like this.

f:=x^2: s:=1/2:
plot(s*f,x=0..4,scaling=constrained,tickmarks=[default,[seq(i=i/s,i=0..8)]]);
f:=x^2+y^3: s:=1/4:
plot3d(s*f,x=0..4,y=0..2,scaling=constrained,tickmarks=[default,default,[seq(i=i/s,i=0..6)]],axes=boxed);
 

Preben Alsholm

f:=Vector(5):
for i to 5 do
f[i] := 5.*i
end do:
f;
 

A:=Matrix(9,9,9);
sqrt~(A); #in Maple 13
map(sqrt,A); #in earlier versions, but also in 13.

B:=Matrix(9,symbol=b);
A/~B; #in Maple 13
zip(`/`,A,B); #in earlier versions, but also in 13.
 

Preben Alsholm

What did you expect to see?

Dirac is a distribution (a generalized function).

Dirac(0) makes no sense and is returned unevaluated by Maple.
Should plot try to evaluate Dirac(t-4) at t = 4 the unevaluated Dirac(0) is just ignored by plot.
This is a nice feature with plot, because it prevents an awful lot of errors in general.

Incidentally, the two and's should be And's, but this would be simpler:

f := piecewise(t < 0, 0,t < 1, 1, t < 2, (-2+t)^2, Dirac(t-4));
 

Preben Alsholm

By using Google I found an implementation of Müller's method by Walter Schreiner at Christian Brothers University:

http://www.cbu.edu/~wschrein/pages/M329Notes.html

Look at Section 2.6: Muller's algorithm, examples are in Muller's method.

Preben Alsholm

add expects finite numbers in the range.

The introduction of k1 is unnecessary it seems to me, and why not define h1 without piecewise and just sum from 0 and up?
restart;
u1:=t->-exp(-a*t);
h1:=t->binomial(n, t)*p^t*(1-p)^(n-t);
Sum(h1(t)*u1(t),t=0..infinity);
value(%);
res:=simplify(%) assuming n::posint;
param:={p=.123,a=.456,n=10};
term:=eval(h1(t)*u1(t),param);
add(term,t=0..1000);
eval(res,param);

However, if you insist

h1:=t->piecewise(t<0,0,binomial(n, t)*p^t*(1-p)^(n-t));

also works with a sum extending from -infinity to infinity.

If you want decimal values (floats), then a, n, and p must have real values. If they are floats, then the result automatically becomes a float. Otherwise use evalf.

Preben Alsholm

Assuming that x0 < x1 < x2 < ...  and that your equations are not equations but just algebraic expressions containing the variable x, you can do:

L:=[seq([x<v[i+1],S[i]],i=0..29)];
p:=piecewise(op(map(op,L))):
plot(p,x=x[0]..x[30]);
 

Preben Alsholm

You forgot a multiplication sign after k. Maple interprets k(xxx) as an application of the function k to xxx.

So do

int((k*(1-x+.3*x*ln(x))+a)^2,x = s .. 1);

Preben Alsholm

GlobalOptimization does not come with Maple 13.

It is an extra you have to buy.

Since I don't have it either, I can reproduce the error message:

with(GlobalOptimization);
Error, invalid input: with expects its 1st argument, pname, to be of type {`module`, package}, but received GlobalOptimization

Preben Alsholm

This may do some of what you want:

restart;
h:=x->sqrt(x);
g:=x->piecewise(type(h(x), And(realcons,Not(infinity) ) ), h(x), 0);
g(-1);
g(2);
g(-infinity);
g(undefined);

##Notice that also a name will return 0:
g(a);
##If that is undesirable you can let the procedure return unevaluated in that case, as in

g:=x->piecewise(type(x,name),'procname(x)', type(h(x),And(realcons,Not(infinity))), h(x), 0);
g(v);
eval(%,v=9);
 

Preben Alsholm

I did this,

restart;
M[Star]:=2.*10^(30);
M[planet]:=2.*10^(27);
q:=(M[planet])/(M[Star]);
#Delta[p]:=max(H,abs(R-a));
Delta4[p]:=max(H^4,(R-a)^4);
Alpha:=unapply(piecewise(R<a, -q^2*6.67*10^(-11)*M[Star]/(2*R)*R^4/Delta4[p] ,R>a,q^2*6.67*10^(-11)*M[Star]/(2*R)*a^4/Delta4[p]) ,R,a);
PDE1:=diff(S(R,t),t)=3/R*diff(R^(1/2)*diff( S(R,t)*R^(1/2)-2*Alpha(R,a)*S(R,t)^(3/2)/(6.67*10^(-11)*M[Star])^(1/2),R),R);
IBC := {S(0, t) = 0, S(99, t) = 0, S(R, 0) = f(R)};
eval(PDE1,{H=1,a=2});
PDE2:=subs(Float(undefined)=0,eval(PDE1,{H=1,a=2}));
#This was done to avoid another error
#The Float(undefined) is caused by differentiating the piecewise defined function.
f:=R->(99-R)*R;
#Just an example of an f
pds1 := pdsolve(PDE2, IBC, numeric, S(R, t), time = t, range = 0 .. 99, spacestep = 0.1e-1, timestep = 0.25e-4);
pds1:-plot(t=0.1,numpoints=5);

 and got the error message

Error, (in pdsolve/numeric/plot) unable to compute solution for t>0.:
unable to store -70000.0000000000+23165.0897338412*I when datatype=float[8]

indicating that S(t) is getting negative (S(t)^(3/2) is in the equation).

Let me add that I have not analyzed your system.

Preben Alsholm
 

You wrote that you want to solve PDE1 for a and H unknown.

Since the pde must be solved numerically, a and H must be given numeric values.

Preben Alsholm

I made these 3 procedures for my own teaching purposes (PolarForm, `print/Exp`, `value/Exp`):

PolarForm:=proc(tal::{algebraic,list,set},e::name:='useexp')
local EXP;
if type(tal,{list,set}) then return map(procname,tal,e) end if;
if type(tal,specfunc(anything,polar))
   then
      if e='useExp' then op(1,tal)*Exp(I*op(2,tal)) else op(1,tal)*'exp'(I*op(2,tal)) end if
   else
procname(radnormal(polar(tal)),e)
end if
end proc:
`print/Exp`:=proc(z) 'e'^(z) end proc:
`value/Exp`:=proc() exp(args) end proc:
PolarForm(1+I*sqrt(3));
%;
PolarForm(1+I*sqrt(3),useExp);
value(%);
 

Preben Alsholm

Maybe the following changed version will help. At least it doesn't produce an error message.

I have used another Delta, Delta4, to avoid abs. Also I have made Alpha a function with unapply

restart;
M[Star]:=2.*10^(30);
M[planet]:=2.*10^(27);
q:=(M[planet])/(M[Star]);
#Delta[p]:=max(H,abs(R-a));
Delta4[p]:=max(H^4,(R-a)^4);
Alpha:=unapply(piecewise(R<a, -q^2*6.67*10^(-11)*M[Star]/(2*R)*R^4/Delta4[p] ,R>a,q^2*6.67*10^(-11)*M[Star]/(2*R)*a^4/Delta4[p]) ,R,a);
PDE1:=diff(S(R,t),t)=3/R*diff(R^(1/2)*diff( S(R,t)*R^(1/2)-2*Alpha(R,a)*S(R,t)^(3/2)/(6.67*10^(-11)*M[Star])^(1/2),R),R);

Preben Alsholm