@wdarrel I'm pretty sure that the numerical algorithms used in solving pdes in Maple do not expect distributions (here Dirac). Whatever comes out is not to be trusted.

Another point: Your pde (this and the original one) can be solved by first making a change of variables like this:
PDEtools:-dchange({x=tau+y,t=tau-y,u(x,t)=w(y,tau)},pde,[y,tau,w]);
#the result is
diff(w(y, tau), tau)+0.591e-1*w(y, tau) = 4.472135955*10^8*Dirac(tau+y-200)+1.0434983895*10^9*Dirac(tau+y-500)
You see that it has the form
(*) diff(w(y, tau), tau)+0.591e-1*w(y, tau) = f(y,tau) 
where f is known. Ignore that f is a distribution and think of it as a function of 2 variables. For every fixed y (*) is a first order ode in tau and has a unique solution if w(y,tau) is known for some tau. Since u(x,0) is known it follows that w(y,tau) is known for tau = y. Thus no more information is needed.

I suggest therefore that you proceed to solve the pde as we solved the original one. Something like this:
restart;
pde := diff(u(x, t), t)+diff(u(x, t), x)+0.591e-1*u(x, t)=447213595.5*Dirac(x-200)+1043498389.5*Dirac(x-500);
ans:=pdsolve(pde);
eval(rhs(ans),t=0)=10^8;
solve(%,{_F1(-x)});
subs(x=x-t,%);
eval(ans,%);
res:=combine(expand(%));
pdetest(res,pde);
eval(res,t=0);
eval(res,x=0);
subs(x=0,res);
U:=collect(expand(rhs(res)),exp(-591/10000*x));
h:=coeff(U,exp(-591/10000*x));
plot3d(h*10^(-21),x=0..1000,t=0..1000,axes=boxed);
plot3d(ln(U),x=0..1000,t=0..1000,axes=boxed);
plot3d(ln(U),x=0..1000,t=0..1,axes=boxed);
plots:-animate(plot,[ln(U),x=0..1000],t=0..20);

@wdarrel You ned to plot the right hand side (rhs) of res and you need plot3d.

Since exponentials are involved the terms get either very large or very small, so logarithmic plots may look more interesting. Also a plot of the coefficient of exp(-x) may be illuminating:

collect(rhs(res),exp);
plot3d(Heaviside(x)-Heaviside(x-t),x=-5..5,t=0..10,axes=boxed);
plot3d(ln(rhs(res)),x=-3..3,t=0..100,axes=boxed);
plot3d(ln(rhs(res)),x=-1..3,t=0..10,axes=boxed);
#An animation
plots:-animate(plot,[ln(rhs(res)),x=-1..3],t=0..10);

@wdarrel You ned to plot the right hand side (rhs) of res and you need plot3d.

Since exponentials are involved the terms get either very large or very small, so logarithmic plots may look more interesting. Also a plot of the coefficient of exp(-x) may be illuminating:

collect(rhs(res),exp);
plot3d(Heaviside(x)-Heaviside(x-t),x=-5..5,t=0..10,axes=boxed);
plot3d(ln(rhs(res)),x=-3..3,t=0..100,axes=boxed);
plot3d(ln(rhs(res)),x=-1..3,t=0..10,axes=boxed);
#An animation
plots:-animate(plot,[ln(rhs(res)),x=-1..3],t=0..10);

Please see the additions above.

Please see the additions above.

You are getting the cpu time. The time it takes the interface to find a way to print the result is not taken into account.

@escorpsy I have uploaded a worksheet containing illustrations and animations. No text, but hopefully somewhat understandable.

MaplePrimes12-11-2.mw

@escorpsy I have uploaded a worksheet containing illustrations and animations. No text, but hopefully somewhat understandable.

MaplePrimes12-11-2.mw

@Markiyan Hirnyk I gave you Brannan and Boyce as a general reference since it wasn't clear to me how much you knew about these things. But you are right, B&B doesn't cover the present situation.

You may want to take a look at

http://www.scholarpedia.org/article/Chetaev_function

After changing variables the present system can be written in vector form as

dx/dt = Dx+|x|e(x), where e(x) ->0 as x->0 and D is a diagonal matrix having the eigenvalues in the diagonal, say in the order -1, 0, 6.

As a Chetaev function you can take V(x1,x2,x3) = -x1^2 + 6*x3^2.

I leave it as an exercise for you to work out the details.

@Markiyan Hirnyk I gave you Brannan and Boyce as a general reference since it wasn't clear to me how much you knew about these things. But you are right, B&B doesn't cover the present situation.

You may want to take a look at

http://www.scholarpedia.org/article/Chetaev_function

After changing variables the present system can be written in vector form as

dx/dt = Dx+|x|e(x), where e(x) ->0 as x->0 and D is a diagonal matrix having the eigenvalues in the diagonal, say in the order -1, 0, 6.

As a Chetaev function you can take V(x1,x2,x3) = -x1^2 + 6*x3^2.

I leave it as an exercise for you to work out the details.

@Markiyan Hirnyk There are quite a few very good textbooks on the subject. The last one I used in an undergraduate class I taught was by Boyce and Brannan.

@Markiyan Hirnyk There are quite a few very good textbooks on the subject. The last one I used in an undergraduate class I taught was by Boyce and Brannan.

I see the point of using ApproximateInt since the integrals are replaced by finite sums and since W(r) and M(r) are just constants in the integration. I give a simple example below.

I notice, that vrad(x) depends on r as well as x, is that intended? Or should (1+r)^(-1/2) have been integrated over too? In your procedure

vproc1 := proc (X, Theta0, Beta, ThetaR) Wsol(parameters = ['theta0' = Theta0, 'beta' = Beta, 'thetaR' = ThetaR]); evalf[15](eval(vradsol, [x = X, 'theta0' = Theta0, 'beta' = Beta, 'thetaR' = ThetaR])) end proc

r doesn't appear.

The use of ApproximateInt seems OK to me. Here is an example.

restart;
eqI := diff(M(r), r) = Int(sqrt(x+M(r)), x = 0 .. M(r));
dsolve({eqI,M(0)=1},numeric); #Doesn't work
diff(eqI,r); #Doesn't seem possible to create a system not involving integrals
#Now ApproximateInt:
eqA := diff(M(r), r) =value( Student:-Calculus1:-ApproximateInt(sqrt(x+M(r)), x = 0 .. M(r),output=sum) );
dsolA:=dsolve({eqA,M(0)=1},numeric,output=listprocedure);
Mn:=subs(dsolA,M(r)):
RI:=subs(M(r)=Mn(r),rhs(eqI));
RA:=subs(M(r)=Mn(r),rhs(eqA)):
plot(RI-RA,r=0..1); #reasonably small

#################

However, to avoid using ApproximateInt you could (certainly in this example) use the procedural form of dsolve/numeric:

solproc := proc(N, r, Y, YP)
    YP[1] := evalf(Int(sqrt(x+Y[1]),x=0..Y[1]));
end proc;
dsol := dsolve(numeric, number=1, procedure=solproc,start=0, initial=Array([1]), procvars=[M(r)]);
plots:-odeplot(dsol,[r,M(r)],0..1);

I see the point of using ApproximateInt since the integrals are replaced by finite sums and since W(r) and M(r) are just constants in the integration. I give a simple example below.

I notice, that vrad(x) depends on r as well as x, is that intended? Or should (1+r)^(-1/2) have been integrated over too? In your procedure

vproc1 := proc (X, Theta0, Beta, ThetaR) Wsol(parameters = ['theta0' = Theta0, 'beta' = Beta, 'thetaR' = ThetaR]); evalf[15](eval(vradsol, [x = X, 'theta0' = Theta0, 'beta' = Beta, 'thetaR' = ThetaR])) end proc

r doesn't appear.

The use of ApproximateInt seems OK to me. Here is an example.

restart;
eqI := diff(M(r), r) = Int(sqrt(x+M(r)), x = 0 .. M(r));
dsolve({eqI,M(0)=1},numeric); #Doesn't work
diff(eqI,r); #Doesn't seem possible to create a system not involving integrals
#Now ApproximateInt:
eqA := diff(M(r), r) =value( Student:-Calculus1:-ApproximateInt(sqrt(x+M(r)), x = 0 .. M(r),output=sum) );
dsolA:=dsolve({eqA,M(0)=1},numeric,output=listprocedure);
Mn:=subs(dsolA,M(r)):
RI:=subs(M(r)=Mn(r),rhs(eqI));
RA:=subs(M(r)=Mn(r),rhs(eqA)):
plot(RI-RA,r=0..1); #reasonably small

#################

However, to avoid using ApproximateInt you could (certainly in this example) use the procedural form of dsolve/numeric:

solproc := proc(N, r, Y, YP)
    YP[1] := evalf(Int(sqrt(x+Y[1]),x=0..Y[1]));
end proc;
dsol := dsolve(numeric, number=1, procedure=solproc,start=0, initial=Array([1]), procvars=[M(r)]);
plots:-odeplot(dsol,[r,M(r)],0..1);

My value for k was chosen so that my dt agreed with yours. Their difference seemed to be zero judging from the plot, simplify didn't do it, but I was satisfied. However, as you will notice, I kept using yours.

I had fun making an animation of the jump and catch: