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-');
.
