Discontinuities when computing integration of error functions using integral function

I am trying to integrate a function over a region in different time intervals. The integration looks something like this.

fun_uz = @(u)1./sqrt(u).*exp(-Z.^2./(2.*u));
fun_Y = @(u)(erf((Y+B)./sqrt(2.*u))-erf((Y-B)./sqrt(2.*u)));
fun_Z = @(u)(erf((X+L+u)./sqrt(2.*u))-erf((X-L+u)./sqrt(2.*u)));
fun = @(u)inc.*fun_uz(u).*fun_Y(u).*fun_Z(u);

fint = integral(fun,0,upperl);

The variable 'upperl' is the upper limit of the integral function. I have to perform this integration over different X,Y, and Z regions and different 'upperl' values. I am getting profiles which are discontinuous for different 'upperl' values. I have shown here profiles at few different 'upperl' values.

I am not able to understand why the discontinuity are occuring, any help is greatly appreciated. Thanks.

[darova]

Attach the whole code

[Yaswanth Sai]

%%

[darova]

I changed line

xs = (-10:1:10)./1000;

Looks ok

[Yaswanth Sai]

Sorry, I forgot to mention, the plots look continuous for the range I have given in the code. They seem to be discontinuous for very large values of 'xs' as shown in the plots I have posted above. For example, for the values mentioned below, the plot is very discontinuous.

xs = (-500:10:10)./1000;
time = 5*10^(-1);

[darova]

time = 1e-2; % The variable which is changed to generate different contour plots
xs = (-500:10:10)./1000;

[Yaswanth Sai]

Yeah, the code works for smaller time values of order 10^(-2) and lesser, but is not working for values greater than that of order 10^(-1) and higher.

Answer 1

[darova]

The function you are trying to integrate looks like following

Put these lines inside for loops

ezplot(fun,[0 upperl])
pause(0.01)

When time > 0.1 upperl is big. When you call integral you don't know how many points it takes. Maybe it misses something

[Yaswanth Sai]

Thank you very much. I understood what you are saying. Is there anything you suggest me to do? I was thinking of subdiving the entire interval until it finds a nonzero value. But this approach does not seem to be computationally efficient.

[Walter Roberson]

Use the waypoints option of integral() if you can.

[Yaswanth Sai]

Yeah, I tried that. The waypoints seem to vary for different values of 'xs'. So I need to figure out a way to assign the limits. Anyways thank you for the suggestion.

[Walter Roberson]

I do not see the whole code posted ?

[Yaswanth Sai]

I added the edited code

Answer 2

[Walter Roberson]

Change the integral to

            fint = integral(fun,0,upperl, 'waypoints', L-X);

You have two erf that only have an input near 0 (and so a measurable output) near-ish -(X+L) to -(X-L) . Some of your integral() calls just happened to evaluate near there, and some of them did not happen to evaluate near there and predicted that there was nothing interesting in that area. The above forces evaluation near that area.

[Yaswanth Sai]

Thank you. I have also used a similar approach to get rid of the discontinuities. I used the following code,

fint = integral(fun,0,upperl,'Waypoints',[abs(X)-abs(X)/2,abs(X)+abs(X)/2]);

Both the approaches are working.