Plotting multivariable numerical functions

Hello. I am new to matlab. I am trying to plot multivariable numeric function. The code is given below:

close all;

%clc;

clear all;

syms k

syms z

syms thetak

syms s1

syms s2

syms s3

syms x

tic

pr=(1250000000000000000000000000000000000*sqrt(2/137))/(4414107*pi^(3/2));zr=pi/1250000;

omega=750000000000000*pi;tau=25*10^(-15);c=3*10^(8);E0=0.0026; %parameters

s21(thetak,k,z,s1,s2,s3) = zr.^2.*tau.*c/omega.*E0.^3.*pi.^(3/2).*1/sqrt(3)...

*exp(-32.*z.^2.*zr.^2/(3.*c.^2.*tau.^2)+1i.*(-4.*s1+2.*s2+2.*s3).*z.*zr.*omega/(3.*c)....

+ 1/3.*1i.*k.*z.*zr.*(1+3.*cos(thetak))-c.*k^.2.*(1+z.^2).*zr.*sin(thetak).^2/(2.*omega.*(3+1i.*(s1-s2-s3).*z)))...

.*((1+1i.*z)/(1-1i.*z)).^(1/2.*(s1-s2-s3))/((3+1i.*(s1-s2-s3).*z).*sqrt(1+z.^2));

s21z(thetak,k,s1,s2,s3) = vpaintegral(s21(thetak,k,z,s1,s2,s3),z,[-Inf,Inf], ...

'ArrayValued',true,'AbsTol',1e-50);

s21zthk(s1,s2,s3) = vpaintegral(cos(thetak/2).^4.*sin(thetak).*vpaintegral(k.^3.*exp(-1/24.*tau.^2.*(c.*k-(s1+s2+s3).*omega).^2)*abs(s21z(thetak,k,s1,s2,s3))^2, ...

k,[1.*10^6, 15.*10^6], ...

'ArrayValued',true,'AbsTol',1e-60),thetak,[0,pi], 'ArrayValued',true);

fplot(s21z(thetak,5.*10^6,1,-1,1),thetak,[0 3.14])

timeElapsed = toc

Here

s21(thetak,k,z,s1,s2,s3) = zr.^2.*tau.*c/omega.*E0.^3.*pi.^(3/2).*1/sqrt(3)...

*exp(-32.*z.^2.*zr.^2/(3.*c.^2.*tau.^2)+1i.*(-4.*s1+2.*s2+2.*s3).*z.*zr.*omega/(3.*c)....

+ 1/3.*1i.*k.*z.*zr.*(1+3.*cos(thetak))-c.*k^.2.*(1+z.^2).*zr.*sin(thetak).^2/(2.*omega.*(3+1i.*(s1-s2-s3).*z)))...

.*((1+1i.*z)/(1-1i.*z)).^(1/2.*(s1-s2-s3))/((3+1i.*(s1-s2-s3).*z).*sqrt(1+z.^2));

is the function I am integrating. and

s21z(thetak,k,s1,s2,s3) = vpaintegral(s21(thetak,k,z,s1,s2,s3),z,[-Inf,Inf], ...

'ArrayValued',true,'AbsTol',1e-50);

is the function where the z integration has been performed and it is a function of the variables thetak, and k. What I would like to do is to plot this function as a function of thetak, for a fixed values of k (the values of s1,s2 and s3 are fixed of course):

fplot(s21z(thetak,5.*10^6,1,-1,1),thetak,[0 3.14])

The code runs and I obtain a blank plot. Can anyone explain why I am getting a blank plot ?

Thank you.

[Torsten]

Do you get a result if you do the integration outside the fplot command ?

E.g.

result = s21z(pi/2,5.*10^6,1,-1,1)

[Taqu3]

Yes I do

double(s21z(pi/2,5*10^6,1,-1,1))

The result I get is:

ans =

1.4479e-48

timeElapsed =

1.0531

[Torsten]

Does the following give a reasonable plot ?

thetak = 0:0.01:3.14;
for i = 1:numel(thetak)
  s21z_array(i) = double(s21z(thetak(i),5*10^6,1,-1,1));
end
plot(thetak,s21z_array)

[Taqu3]

Yes and simple enough that I can understand :) Thank you!