% This script calculates the resonances for square-well potential with % additional barrier % Values

% This script calculates the resonances for square-well potential with

% additional barrier

% Values clear all

[maxit,tol,R,Rmax,V0]=deal(200,1e-8,0.01:0.01:0.5,0.5,300);

% Functions syms E R1 m;

y(m,E,R1) = besselj(m,sqrt(E))*(sqrt(E-V0)*bessely(m,sqrt(E)*R1)*besselj(m-1,sqrt(E-V0)*R1)-sqrt (E)*besselj(m,sqrt(E-V0)*R1)*bessely(m-1,sqrt(E)*R1))-bessely(m,sqrt(E))*(sqrt(E-V0)*besselj(m,sqrt (E)*R1)*besselj(m-1,sqrt(E-V0)*R1)-sqrt(E)*besselj(m,sqrt(E-V0)*R1)*besselj(m-1,sqrt(E)*R1));

yprime(m,E,R1) = diff(y,E);

f = matlabFunction(y);

fprime = matlabFunction(yprime);

% Initial guess matrix

% This uses the besselzero script avaiable at http://nl.mathworks.com/matlabcentral/fileexchange/

submissions/6794/v/1/download/zip x=besselzero(0,40,1);

for m=1:1:20

x=[x,besselzero(m,40,1)];

end

% Perturb for better results x=x.^2-0.01-i/100;

disp('Made matrix, now solving.') for m=0:1:20 % Loop over m

for r=R % increase the new barrier for j=1:length(x) % for each zero

for i=1:maxit% iterate for some zero xold=x(j,m+1);

x(j,m+1)=x(j,m+1)-f(m,x(j,m+1),r)/fprime(m,x(j,m+1),r);

err=xold-x(j,m+1);

if abs(err)<tol break end; end; end; end; end

% Testing the found values test=1;

if test

disp('Now testing the values.') for m=0:1:20

for j=1:40

a=f(m,x(j,m+1),Rmax);

if abs(a)>1e-12;

warning('Too high') [a,x(j,m+1)]

break end; end; end; end

% Plotting the values plot=1;

if plot

disp('Making a plot.') figure

for m=0:1:20

y=real(x(:,m+1));

m_vec=m*ones(1,40);

scatter(m_vec,y,'filled');set(gca,'FontSize',15) hold on

scatter(-m_vec,y,'filled') end

axis([-20 20 0 500]) xlabel('m')

ylabel('E')

str=sprintf('Spectrum for R1=%f',Rmax);

title(str) scatter(0,V0) end