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
 

Say your vector is

v:= <a(t),b(t),c(t)>;

then in Maple 13 you can differentiate elementwise like this

diff~(v. t):

where you should notice the tilde.

In earlier versions:

map(diff, v, t);

map is still useful, though.

Preben Alsholm

Your procedure with two changes, viz. the local declaration and contourplot replaced by its long name:

PlotMfield := proc (a::float, b::float, c::float, d::float, p::float, chi::float)
local g,G,f,Bz,Bphi,Br;
g := z -> c*(exp(-(1/2)*z^2/a^2)/a-exp(-(1/2)*z^2/b^2)/b)/sqrt(2*Pi) ;
G :=z -> int(g(zeta), zeta = -infinity .. z);
f  := r -> d*exp(-p*r);
Bz := chi*p*(p-1/r)*G(z)*f(r);
Bphi := p*g(z)*f(r);
Br := chi*p*g(z)*f(r);
plots:-contourplot(sqrt(Br^2+Bphi^2+Bz^2), r = 0 .. 15, z = (-1)*2.5 .. 2.5, grid = [15, 75], scaling = constrained);
end proc;
 

Your second problem can be solved by executing
printlevel:=2; 
(because of the nested loop)
or by giving an explicit print command as in

for i in La do
   for j in Lc do
     print(PlotMfield (i,1.1*i,j,1.0,0.4,1.0))
   end do;
end do;

Preben Alsholm
 

The integral is never zero, and you need to change 'pi' to 'Pi'.

You can try the following instead.

F:=(t,E)->Int(sqrt((1-(1-E^2)*sin(x)^2)/(1-(1-E^2)*t^2*sin(x)^2)),x=0..Pi/2);
plots:-contourplot(F(t,E),t=0..1,E=0..1);
 kk:=evalf(F(.3,.4));
plots:-implicitplot(F(t,E)=kk,t=0..1,E=0..1);
 

Preben Alsholm

You could use concatenation, i.e. cat in Maple, as in

cat("Figure ", j+1);

I couldn't test it on your procedure since several things are missing (and you seem to be solving for a constant in b:=[solve(t^3-a=0,a)]);

Preben Alsholm