Wrong plot upon right calculation

Hi,

I have the next problem. Here is a code for calculation of the function K(s) with integral fun1. When s=0 the K(s)=1 at all values of parameters n,t,r1. And when I calculate it by vpa for s=0 I obtain this K=1 at all parameters. Thus, the curve K(s) must always start from the point (0,1). But when I make a plot, the curve start from some different value. Please, help me: what do I do wrong? Calculation of this integral by syms package gives right plots, but it calculates too long and I want to do it by the code below. Thank you.

clear all, close all

D=1;

n=0.01;

r1=1;

t=1;

r=(r1^2+4*t/3);

%s=0;

s=0:0.1:10;

for i = 1:length(s)

k=s(i);

fun1 = @(x)((((1-(((1-x.^2).*k.^2)/(6*r))).*besseli(0,(((1-x.^2).*k.^2)/(6*r))/2))+(((((1-x.^2).*k.^2)/(6*r))).*besseli(1,(((1-x.^2).*k.^2)/(6*r))/2))).*exp(((k.^2)/(12*r)-2*n*t).*x.^2));

f1(i,:)=integral(fun1,-1,1);

K=(exp(-k.^2/(12*r))*D*sqrt(2*n*t)/(sqrt(pi)*erf(sqrt(2*n*t))))*f1;

end

plot(s,K,'b-');

%vpa(Fun,5) % calculation of K(s,t) at certain s. At s=0 must be next K(s,t)=1

Answer 1

[Star Strider]

If it needs to begin at (0,1), the only changes required are that in the ‘K’ assignment, the current value of ‘f1’ needs to be referenced specifically and ‘K’ needs to be indexed:

K(i) = (exp(-k.^2/(12*r))*D*sqrt(2*n*t)/(sqrt(pi)*erf(sqrt(2*n*t))))*f1(i);

See if those changes produce the result you want —

D=1;
n=0.01;
r1=1;
t=1;
r=(r1^2+4*t/3);
%s=0;
s=0:0.1:10;
for i = 1:length(s)   
    k=s(i);
    fun1 = @(x)((((1-(((1-x.^2).*k.^2)/(6*r))).*besseli(0,(((1-x.^2).*k.^2)/(6*r))/2))+(((((1-x.^2).*k.^2)/(6*r))).*besseli(1,(((1-x.^2).*k.^2)/(6*r))/2))).*exp(((k.^2)/(12*r)-2*n*t).*x.^2));
    f1(i,:)=integral(fun1,-1,1);  
    K(i)=(exp(-k.^2/(12*r))*D*sqrt(2*n*t)/(sqrt(pi)*erf(sqrt(2*n*t))))*f1(i);      
end
plot(s,K,'b-');

.

[Hexe]

Dear Star Strider,

Thank you very much for you help. Now it works as it has to.

Sincerely

Olha.

[Star Strider]

As always, my pleasure!