@acer Hi! When I check your answer again, I realized I could not understand the following part:

alt_c3_1 := simplify(combine(student[changevar](Y=1/(y-1), c3_1, Y)),size) assuming r>0, u>0, u<10:
alt_c3_1 := Int(op(1,alt_c3_1), Y=1..100, epsilon=1e-6, method=_d01akc);

the change of variable is Y=1/(y-1) and integral with y dummy variable is from 1 to 2. With new Y dummy variable, integral should be from infinity (can be 100) to 1. But you considered   Y=1..100. Actually my question is what does op(1,alt_c3_1) mean? I would like to apply this trick for my other integrals and I want to understand it.

Thank you!

@acer  Thank you for your suggestion about method=ftocms. In physical bacground, the figure should be constant. How can I use the medhof which you suggested?

@Axel Vogt  I checked the equation, everything is correct. Because of convolution theorem, we can change only the argument order as follow:

      c2:=f1(y)*(sin((x-y)^2/8/r+Pi/4)+sin((x+y)^2/8/r+Pi/4))*((((x+y)-4*d1*(t1-r))/(4*Pi*(t1-r)*sqrt((t1-r)*r)))*(exp(-(x+y)^2/(8*(t1-r))))+((2*d1^2*sqrt(1/(2*Pi*r)))*exp(d1*(x+y)+(2*d1^2*(t1-r))))*erfc(((x+y)+4*d1*(t1-r))/(sqrt(8*(t1-r)))));

@tomleslie Thank you! I did not expect they were very close to zero. I will check my equations..

Edit: Hi Tomleslie! I checked my equations and there are correct. I do not know why the results are so close to zero. Because I need to change the value of d1 as 0.2, 0.4, 1 and 2 in order to see the effectiveness of d1. Only I can change the order of arguments inside of the integral because of the convolution theorem.

c2:=f1(y)*(sin((x-y)^2/8/r+Pi/4)+sin((x+y)^2/8/r+Pi/4))*((((x+y)-4*d1*(t1-r))/(4*Pi*(t1-r)*sqrt((t1-r)*r)))*(exp(-(x+y)^2/(8*(t1-r))))+((2*d1^2*sqrt(1/(2*Pi*r)))*exp(d1*(x+y)+(2*d1^2*(t1-r))))*erfc(((x+y)+4*d1*(t1-r))/(sqrt(8*(t1-r)))));

When I check the points, it gives

u2_pre := Int(c2, r= 0..0.2, numeric):
u21_pre := Int(u2_pre, y=0..100, numeric):
u21:= unapply(u21_pre,x):
seq( value(u21(j)), j=0..10);

 Float(undefined), Float(undefined), Float(undefined),

                  -15                -16                 -20  
   -1.716493742 10   , 2.915266980 10   , -2.970490519 10   ,

                 -25                -30                -36  
   8.672301117 10   , 7.964105024 10   , 2.255831187 10   ,

                  -43                 -49
   -8.859604626 10   , -4.786537728 10   

How can I get rid of the Floating problem?

@tomleslie The problem is that it should not be zero.

@Carl Love : Thank you it helped alot. Now I need to understand why I do not have constant figure..Thanks again!

@Carl Love : Thank you for your answer! Interesting approach.. As I understand I can use G50 as g5. But when I substitute those codes , I cannot plot E2..

Thank you guys for your answers! I could plot g1 with using heaviside f1 function.E2 should give constant figure because I would like to show energy is conserved, that is no energy changing. If I change f1 with

f1:=y->exp(-(y-4)^2/1):

Then

E2(0);
Error, (in E2) numeric exception: division by zero.

I do not know how I can plot it.

@acer : Thank you for your answer! I aggre with you about smoothness, but f1(y) should be a heaviside function. I just improved the mw.file. y interval is from 0 to infinity. But for numerical integration I cannot consider it bigger than 42 for u2, u21 and bigger than 10 for u22 and u23 . I can plot for t=0, t=0.2 and t=0.4, but I need t=1 and t=2 as well. Somehow it blows up and I cannot plot them together. do you know how I can?

P.S. I think for increasing time I need to divide y interval for smoothness.

tryingf.mw

@acer  Hi, I liked the way which you tried to solve C33. You used change of variable in (1,2) interval.

I would like to use the same trick to the new C33 in (0,t1) interval as follow, but I could not understand why the integral gives zero for every value?

restart: plots:-setoptions(gridlines=false);
H := a -> piecewise(a>=0,1):
f1 := x->(H(x-1*Pi)-H(x-2*Pi))*sin(x)^2:
d1:=1:t1:=2/10:
c33:=f1(y)*(sin((x-y)^2/8/(t1-r)+Pi/4)+sin((x+y)^2/8/(t1-r)+Pi/4))*((((x+y)-4*d1*r)/(4*Pi*r*sqrt(r*(t1-r))))*(exp(-(x+y)^2/(8*t1)))+((2*d1^2*sqrt(1/(2*Pi*(t1-r))))*exp(d1*(x+y)+(2*d1^2*r))*erfc(((x+y)+4*d1*r)/(2*sqrt(2*r))))):

I enclosed the file tryingf.mw

If you can help me I would be happy!

Thanks

@acer  I am so happy , thank you for your answer! I learn a lot from the answer. Now I am trying to understand all progress. But, I had a quick look. It seems perfect for me! Well done! 

@Preben Alsholm Thank you so much for showing the typo. I have just edited it again..