How can I plot the following function?

alpha=0.001:0.1:5;
s=1;sb=.1;s1=sb./s;j=0.1;EZ=1;zeta0=-1;beta=1;i=(-1).^(1./2);
sigma=0.1; ur=1;k=5;b1=0.1;b2=0.1;c=0.1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
k1=k.^2.*zeta0./(1+c);
T1=(1+c)./(2.*s.*(i.*alpha.^2.*j-4.*c));
T2=(4.*s.*c-i.*alpha.^2)./(1+c);
r11=(4.*s.*c+i.*alpha.^2.*(c-j.*s.*(1+c)))./(c.*(1+c));r22=alpha.^2.*s.*(4.*i.*c+alpha.^2.*j)./(c.*(1+c));
alpha11=((1./2).*(r11+(r11.^2-4.*r22).^(1./2))).^(1./2);alpha12=-((1./2).*(r11+(r11.^2-4.*r22).^(1./2))).^(1./2);
alpha22=((1./2).*(r11-(r11.^2-4.*r22).^(1./2))).^(1./2);alpha21=-((1./2).*(r11-(r11.^2-4.*r22).^(1./2))).^(1./2);
KK=-(k1.*(k.^2-s.*(4-i.*alpha.^2.*j.*c.^(-1))))./(k.^4-(r11).*k.^2+r22);
format long g
syms r
L3 = vpaintegral(r*besselk(1,alpha11*r),r,sigma,1);
w1 =2.*( -T1 .* alpha11 .* L3.*(alpha11 .^ 2 + T2));
plot(alpha,abs(w1),'-black','LineWidth',1.2);
ylabel('$w1(r)$','Interpreter','latex','FontSize',12,'FontName','TiemsNewRoman','FontWeight','Normal')
xlabel('$\alpha$','Interpreter','latex','FontSize',12,'FontName','TiemsNewRoman','FontWeight','Normal')

%%%%%%%%%%%

[Torsten]

I get a plot, but note that the elements of L3 and thus of w1 are complex. Thus you can plot abs(w1), real(w1) or imag(w1), but not w1 itself.

[Shreen El-Sapa]

Yes, no problem

[Torsten]

So problem solved ?

[Shreen El-Sapa]

where? can you excute it? I couldnot get any thing.

[Torsten]

See above. I reduced the alpha resolution from 0.0001 to 0.1 to save computation time (and I think this resolution is sufficient).

[Shreen El-Sapa]

Thanks so much.