% Computes the largest eigenvalue in the para case.
% To create a table, type, e.g.:
% >>  for Z=1:10, for om=8:11, tab(Z,om-7) = PekerisPara(om,Z); end; end
% To see the result, type:
% >>  format long; tab
% To write this table to a file in the format as it appears in the paper:
% >>  dlmwrite('C:\test.txt', tab, 'delimiter', ' ', 'precision', '%0.10f');

function ev = PekerisPara(omega, Z)
dim = ceil((omega+1)*(omega+3)*(2*omega+7)/24);
rowi = ones(1,33*dim); 
coli = ones(1,33*dim); 
valsA = zeros(1,33*dim);
valsB = zeros(1,33*dim);
i = 1;

for w=0:omega, 
    for lm=0:w, 
        for l=0:floor(lm/2),
            m = lm-l;
            n = w-lm;
            kr = (2*w*(w+2)*(2*w+5) + 3*(1-(-1)^w) + 12*(l+m)^2 + 6*(1-(-1)^(l+m)) + 24*(l+m))/48 + l+1;

			% term no. 1
			if l>=2,
				if l<=2+m,
					kc = (-39-6*(-1)^(l+m)-3*(-1)^w+12*l^2-24*m+12*m^2+24*l*(1+m)-4*w-6*w^2+4*w^3)/48;
				else
					kc = (57-6*(-1)^(l+m)-3*(-1)^w+12*l^2+24*l*(-1+m)+24*m+12*m^2-4*w-6*w^2+4*w^3)/48;
				end
				if kc <= dim,
					rowi(i) = kr; coli(i) = kc;
					valsA(i) = 4*(-1+l)*l*Z;
					valsB(i) = 4*(-1+l)*l*(1+m+n);
					i = i+1;
				end
			end

			% term no. 2
			if l>=2,
				if l<=2+m,
					kc = (-45-6*(-1)^(l+m)+3*(-1)^w+12*l^2-24*m+12*m^2+24*l*(1+m)-4*w+6*w^2+4*w^3)/48;
				else
					kc = (51-6*(-1)^(l+m)+3*(-1)^w+12*l^2+24*l*(-1+m)+24*m+12*m^2-4*w+6*w^2+4*w^3)/48;
				end
				if kc <= dim,
					rowi(i) = kr; coli(i) = kc;
					valsA(i) = 0;
					valsB(i) = 4*(-1+l)*l*(1+n);
					i = i+1;
				end
			end

			% term no. 3
			if l>=1&&m>=1,
				if l<=m,
					kc = (9-6*(-1)^(l+m)-3*(-1)^w+12*l^2-24*m+12*m^2+24*l*(1+m)-4*w-6*w^2+4*w^3)/48;
				else
					kc = (9-6*(-1)^(l+m)-3*(-1)^w+12*l^2+24*l*(-1+m)+24*m+12*m^2-4*w-6*w^2+4*w^3)/48;
				end
				if kc <= dim,
					rowi(i) = kr; coli(i) = kc;
					valsA(i) = 4*l*m*(-1+2*Z);
					valsB(i) = 4*l*m*(l+m);
					i = i+1;
				end
			end

			% term no. 4
			if l>=1&&n>=1,
				if l<=1+m,
					kc = (-3+6*(-1)^(l+m)-3*(-1)^w+12*l^2+12*m^2+24*l*(2+m)-4*w-6*w^2+4*w^3)/48;
				else
					kc = (45+6*(-1)^(l+m)-3*(-1)^w+12*l^2+48*m+24*l*m+12*m^2-4*w-6*w^2+4*w^3)/48;
				end
				if kc <= dim,
					rowi(i) = kr; coli(i) = kc;
					valsA(i) = 2*l*n*(-1+2*Z);
					valsB(i) = 2*l*n*(1+2*m+n);
					i = i+1;
				end
			end

			% term no. 5
			if l>=1,
				if l<=1+m,
					kc = (-9+6*(-1)^(l+m)+3*(-1)^w+12*l^2+12*m^2+24*l*(2+m)-4*w+6*w^2+4*w^3)/48;
				else
					kc = (39+6*(-1)^(l+m)+3*(-1)^w+12*l^2+48*m+24*l*m+12*m^2-4*w+6*w^2+4*w^3)/48;
				end
				if kc <= dim,
					rowi(i) = kr; coli(i) = kc;
					valsA(i) = -2*l*(-3+6*Z+8*l*Z+n*(-2+4*Z)+m*(-4+8*Z));
					valsB(i) = -2*l*(2+4*m^2+3*n+n^2+12*l*(1+m+n)+m*(6+4*n));
					i = i+1;
				end
			end

			% term no. 6
			if l>=1,
				if l<=1+m,
					kc = (-3+6*(-1)^(l+m)-3*(-1)^w+12*l^2+12*m^2+24*l*(2+m)+20*w+18*w^2+4*w^3)/48;
				else
					kc = (45+6*(-1)^(l+m)-3*(-1)^w+12*l^2+48*m+24*l*m+12*m^2+20*w+18*w^2+4*w^3)/48;
				end
				if kc <= dim,
					rowi(i) = kr; coli(i) = kc;
					valsA(i) = 2*l*(1+n)*(-1+2*Z);
					valsB(i) = -2*l*(1+n)*(4*l-2*m+n);
					i = i+1;
				end
			end

			% term no. 7
			if l>=1,
				if l<=1+m,
					kc = (39+6*(-1)^(l+m)+3*(-1)^w+12*l^2+12*m^2+24*l*(2+m)+68*w+30*w^2+4*w^3)/48;
				else
					kc = (87+6*(-1)^(l+m)+3*(-1)^w+12*l^2+48*m+24*l*m+12*m^2+68*w+30*w^2+4*w^3)/48;
				end
				if kc <= dim,
					rowi(i) = kr; coli(i) = kc;
					valsA(i) = 0;
					valsB(i) = 2*l*(2+3*n+n^2);
					i = i+1;
				end
			end

			% term no. 8
			if l>=1,
				if l<=2+m,
					kc = (9-6*(-1)^(l+m)-3*(-1)^w+12*l^2+24*m+12*m^2+24*l*(3+m)+20*w+18*w^2+4*w^3)/48;
				else
					kc = (105-6*(-1)^(l+m)-3*(-1)^w+12*l^2+72*m+12*m^2+24*l*(1+m)+20*w+18*w^2+4*w^3)/48;
				end
				if kc <= dim,
					rowi(i) = kr; coli(i) = kc;
					valsA(i) = 4*l*(1+m)*(-1+2*Z);
					valsB(i) = 4*l*(1+m)*(1+l+m);
					i = i+1;
				end
			end

			% term no. 9
			if m>=2,
				if 2+l<=m,
					kc = (57-6*(-1)^(l+m)-3*(-1)^w+12*l^2-24*m+12*m^2+24*l*(1+m)-4*w-6*w^2+4*w^3)/48;
				else
					kc = (-39-6*(-1)^(l+m)-3*(-1)^w+12*l^2+24*l*(-1+m)+24*m+12*m^2-4*w-6*w^2+4*w^3)/48;
				end
				if kc <= dim,
					rowi(i) = kr; coli(i) = kc;
					valsA(i) = 4*(-1+m)*m*Z;
					valsB(i) = 4*(-1+m)*m*(1+l+n);
					i = i+1;
				end
			end

			% term no. 10
			if m>=2,
				if 2+l<=m,
					kc = (51-6*(-1)^(l+m)+3*(-1)^w+12*l^2-24*m+12*m^2+24*l*(1+m)-4*w+6*w^2+4*w^3)/48;
				else
					kc = (-45-6*(-1)^(l+m)+3*(-1)^w+12*l^2+24*l*(-1+m)+24*m+12*m^2-4*w+6*w^2+4*w^3)/48;
				end
				if kc <= dim,
					rowi(i) = kr; coli(i) = kc;
					valsA(i) = 0;
					valsB(i) = 4*(-1+m)*m*(1+n);
					i = i+1;
				end
			end

			% term no. 11
			if m>=1&&n>=1,
				if 1+l<=m,
					kc = (45+6*(-1)^(l+m)-3*(-1)^w+12*l^2+12*m^2+24*l*(2+m)-4*w-6*w^2+4*w^3)/48;
				else
					kc = (-3+6*(-1)^(l+m)-3*(-1)^w+12*l^2+48*m+24*l*m+12*m^2-4*w-6*w^2+4*w^3)/48;
				end
				if kc <= dim,
					rowi(i) = kr; coli(i) = kc;
					valsA(i) = 2*m*n*(-1+2*Z);
					valsB(i) = 2*m*n*(1+2*l+n);
					i = i+1;
				end
			end

			% term no. 12
			if m>=1,
				if 1+l<=m,
					kc = (39+6*(-1)^(l+m)+3*(-1)^w+12*l^2+12*m^2+24*l*(2+m)-4*w+6*w^2+4*w^3)/48;
				else
					kc = (-9+6*(-1)^(l+m)+3*(-1)^w+12*l^2+48*m+24*l*m+12*m^2-4*w+6*w^2+4*w^3)/48;
				end
				if kc <= dim,
					rowi(i) = kr; coli(i) = kc;
					valsA(i) = -2*m*(-3+6*Z+8*m*Z+n*(-2+4*Z)+l*(-4+8*Z));
					valsB(i) = -2*m*(4*l^2+(1+n)*(2+12*m+n)+2*l*(3+6*m+2*n));
					i = i+1;
				end
			end

			% term no. 13
			if m>=1,
				if 1+l<=m,
					kc = (45+6*(-1)^(l+m)-3*(-1)^w+12*l^2+12*m^2+24*l*(2+m)+20*w+18*w^2+4*w^3)/48;
				else
					kc = (-3+6*(-1)^(l+m)-3*(-1)^w+12*l^2+48*m+24*l*m+12*m^2+20*w+18*w^2+4*w^3)/48;
				end
				if kc <= dim,
					rowi(i) = kr; coli(i) = kc;
					valsA(i) = 2*m*(1+n)*(-1+2*Z);
					valsB(i) = -2*m*(1+n)*(-2*l+4*m+n);
					i = i+1;
				end
			end

			% term no. 14
			if m>=1,
				if 1+l<=m,
					kc = (87+6*(-1)^(l+m)+3*(-1)^w+12*l^2+12*m^2+24*l*(2+m)+68*w+30*w^2+4*w^3)/48;
				else
					kc = (39+6*(-1)^(l+m)+3*(-1)^w+12*l^2+48*m+24*l*m+12*m^2+68*w+30*w^2+4*w^3)/48;
				end
				if kc <= dim,
					rowi(i) = kr; coli(i) = kc;
					valsA(i) = 0;
					valsB(i) = 2*m*(2+3*n+n^2);
					i = i+1;
				end
			end

			% term no. 15
			if n>=2,
				if l<=m,
					kc = (57-6*(-1)^(l+m)-3*(-1)^w+12*l^2+24*m+12*m^2+24*l*(3+m)-4*w-6*w^2+4*w^3)/48;
				else
					kc = (57-6*(-1)^(l+m)-3*(-1)^w+12*l^2+72*m+12*m^2+24*l*(1+m)-4*w-6*w^2+4*w^3)/48;
				end
				if kc <= dim,
					rowi(i) = kr; coli(i) = kc;
					valsA(i) = -((-1+n)*n);
					valsB(i) = 0;
					i = i+1;
				end
			end

			% term no. 16
			if n>=1,
				if l<=m,
					kc = (51-6*(-1)^(l+m)+3*(-1)^w+12*l^2+24*m+12*m^2+24*l*(3+m)-4*w+6*w^2+4*w^3)/48;
				else
					kc = (51-6*(-1)^(l+m)+3*(-1)^w+12*l^2+72*m+12*m^2+24*l*(1+m)-4*w+6*w^2+4*w^3)/48;
				end
				if kc <= dim,
					rowi(i) = kr; coli(i) = kc;
					valsA(i) = 4*n*(1+l+m+n-2*Z-2*l*Z-2*m*Z);
					valsB(i) = -4*n*(l-l^2+m+4*l*m-m^2+2*n+2*l*n+2*m*n);
					i = i+1;
				end
			end

			% term no. 17
			if l<=m,
				kc = (57-6*(-1)^(l+m)-3*(-1)^w+12*l^2+24*m+12*m^2+24*l*(3+m)+20*w+18*w^2+4*w^3)/48;
			else
				kc = (57-6*(-1)^(l+m)-3*(-1)^w+12*l^2+72*m+12*m^2+24*l*(1+m)+20*w+18*w^2+4*w^3)/48;
			end
			if kc <= dim,
				rowi(i) = kr; coli(i) = kc;
				valsA(i) = 2*(-5-7*n-3*n^2+16*Z+12*l^2*Z+12*m^2*Z+8*n*Z+m*(-6-4*n+24*Z+8*n*Z)+2*l*(-3-2*n+12*Z+4*n*Z+m*(-4+8*Z)));
				valsB(i) = 8*(5*m^2*(1+n)+5*l^2*(1+m+n)+2*(2+3*n+n^2)+m*(8+9*n+2*n^2)+l*(8+5*m^2+9*n+2*n^2+4*m*(3+n)));
				i = i+1;
			end

			% term no. 18
			if l<=m,
				kc = (99-6*(-1)^(l+m)+3*(-1)^w+12*l^2+24*m+12*m^2+24*l*(3+m)+68*w+30*w^2+4*w^3)/48;
			else
				kc = (99-6*(-1)^(l+m)+3*(-1)^w+12*l^2+72*m+12*m^2+24*l*(1+m)+68*w+30*w^2+4*w^3)/48;
			end
			if kc <= dim,
				rowi(i) = kr; coli(i) = kc;
				valsA(i) = 4*(1+n)*(2+l+m+n-2*Z-2*l*Z-2*m*Z);
				valsB(i) = 4*(1+n)*(l^2+m^2-2*(1+n)-m*(3+2*n)-l*(3+4*m+2*n));
				i = i+1;
			end

			% term no. 19
			if l<=m,
				kc = (201-6*(-1)^(l+m)-3*(-1)^w+12*l^2+24*m+12*m^2+24*l*(3+m)+140*w+42*w^2+4*w^3)/48;
			else
				kc = (201-6*(-1)^(l+m)-3*(-1)^w+12*l^2+72*m+12*m^2+24*l*(1+m)+140*w+42*w^2+4*w^3)/48;
			end
			if kc <= dim,
				rowi(i) = kr; coli(i) = kc;
				valsA(i) = -2-3*n-n^2;
				valsB(i) = 0;
				i = i+1;
			end

			% term no. 20
			if n>=2,
				if l<=1+m,
					kc = (87+6*(-1)^(l+m)+3*(-1)^w+12*l^2+48*m+12*m^2+24*l*(4+m)-4*w+6*w^2+4*w^3)/48;
				else
					kc = (135+6*(-1)^(l+m)+3*(-1)^w+12*l^2+96*m+12*m^2+24*l*(2+m)-4*w+6*w^2+4*w^3)/48;
				end
				if kc <= dim,
					rowi(i) = kr; coli(i) = kc;
					valsA(i) = 0;
					valsB(i) = 2*(1+m)*(-1+n)*n;
					i = i+1;
				end
			end

			% term no. 21
			if n>=1,
				if l<=1+m,
					kc = (93+6*(-1)^(l+m)-3*(-1)^w+12*l^2+48*m+12*m^2+24*l*(4+m)+20*w+18*w^2+4*w^3)/48;
				else
					kc = (141+6*(-1)^(l+m)-3*(-1)^w+12*l^2+96*m+12*m^2+24*l*(2+m)+20*w+18*w^2+4*w^3)/48;
				end
				if kc <= dim,
					rowi(i) = kr; coli(i) = kc;
					valsA(i) = 2*(1+m)*n*(-1+2*Z);
					valsB(i) = -2*(1+m)*n*(3-2*l+4*m+n);
					i = i+1;
				end
			end

			% term no. 22
			if l<=1+m,
				kc = (135+6*(-1)^(l+m)+3*(-1)^w+12*l^2+48*m+12*m^2+24*l*(4+m)+68*w+30*w^2+4*w^3)/48;
			else
				kc = (183+6*(-1)^(l+m)+3*(-1)^w+12*l^2+96*m+12*m^2+24*l*(2+m)+68*w+30*w^2+4*w^3)/48;
			end
			if kc <= dim,
				rowi(i) = kr; coli(i) = kc;
				valsA(i) = -2*(1+m)*(-3+14*Z+8*m*Z+n*(-2+4*Z)+l*(-4+8*Z));
				valsB(i) = -2*(1+m)*(4*l^2+(1+n)*(14+12*m+n)+2*l*(9+6*m+2*n));
				i = i+1;
			end

			% term no. 23
			if l<=1+m,
				kc = (237+6*(-1)^(l+m)-3*(-1)^w+12*l^2+48*m+12*m^2+24*l*(4+m)+140*w+42*w^2+4*w^3)/48;
			else
				kc = (285+6*(-1)^(l+m)-3*(-1)^w+12*l^2+96*m+12*m^2+24*l*(2+m)+140*w+42*w^2+4*w^3)/48;
			end
			if kc <= dim,
				rowi(i) = kr; coli(i) = kc;
				valsA(i) = 2*(1+m)*(1+n)*(-1+2*Z);
				valsB(i) = 2*(1+m)*(1+n)*(2+2*l+n);
				i = i+1;
			end

			% term no. 24
			if n>=1,
				if l<=2+m,
					kc = (195-6*(-1)^(l+m)+3*(-1)^w+12*l^2+72*m+12*m^2+24*l*(5+m)+68*w+30*w^2+4*w^3)/48;
				else
					kc = (291-6*(-1)^(l+m)+3*(-1)^w+12*l^2+120*m+12*m^2+24*l*(3+m)+68*w+30*w^2+4*w^3)/48;
				end
				if kc <= dim,
					rowi(i) = kr; coli(i) = kc;
					valsA(i) = 0;
					valsB(i) = 4*(2+3*m+m^2)*n;
					i = i+1;
				end
			end

			% term no. 25
			if l<=2+m,
				kc = (297-6*(-1)^(l+m)-3*(-1)^w+12*l^2+72*m+12*m^2+24*l*(5+m)+140*w+42*w^2+4*w^3)/48;
			else
				kc = (393-6*(-1)^(l+m)-3*(-1)^w+12*l^2+120*m+12*m^2+24*l*(3+m)+140*w+42*w^2+4*w^3)/48;
			end
			if kc <= dim,
				rowi(i) = kr; coli(i) = kc;
				valsA(i) = 4*(2+3*m+m^2)*Z;
				valsB(i) = 4*(2+3*m+m^2)*(1+l+n);
				i = i+1;
			end

			% term no. 26
			if m>=1,
				if 2+l<=m,
					kc = (105-6*(-1)^(l+m)-3*(-1)^w+12*l^2+24*m+12*m^2+24*l*(3+m)+20*w+18*w^2+4*w^3)/48;
				else
					kc = (9-6*(-1)^(l+m)-3*(-1)^w+12*l^2+72*m+12*m^2+24*l*(1+m)+20*w+18*w^2+4*w^3)/48;
				end
				if kc <= dim,
					rowi(i) = kr; coli(i) = kc;
					valsA(i) = 4*(1+l)*m*(-1+2*Z);
					valsB(i) = 4*(1+l)*m*(1+l+m);
					i = i+1;
				end
			end

			% term no. 27
			if n>=2,
				if 1+l<=m,
					kc = (135+6*(-1)^(l+m)+3*(-1)^w+12*l^2+48*m+12*m^2+24*l*(4+m)-4*w+6*w^2+4*w^3)/48;
				else
					kc = (87+6*(-1)^(l+m)+3*(-1)^w+12*l^2+96*m+12*m^2+24*l*(2+m)-4*w+6*w^2+4*w^3)/48;
				end
				if kc <= dim,
					rowi(i) = kr; coli(i) = kc;
					valsA(i) = 0;
					valsB(i) = 2*(1+l)*(-1+n)*n;
					i = i+1;
				end
			end

			% term no. 28
			if n>=1,
				if 1+l<=m,
					kc = (141+6*(-1)^(l+m)-3*(-1)^w+12*l^2+48*m+12*m^2+24*l*(4+m)+20*w+18*w^2+4*w^3)/48;
				else
					kc = (93+6*(-1)^(l+m)-3*(-1)^w+12*l^2+96*m+12*m^2+24*l*(2+m)+20*w+18*w^2+4*w^3)/48;
				end
				if kc <= dim,
					rowi(i) = kr; coli(i) = kc;
					valsA(i) = 2*(1+l)*n*(-1+2*Z);
					valsB(i) = -2*(1+l)*n*(3+4*l-2*m+n);
					i = i+1;
				end
			end

			% term no. 29
			if 1+l<=m,
				kc = (183+6*(-1)^(l+m)+3*(-1)^w+12*l^2+48*m+12*m^2+24*l*(4+m)+68*w+30*w^2+4*w^3)/48;
			else
				kc = (135+6*(-1)^(l+m)+3*(-1)^w+12*l^2+96*m+12*m^2+24*l*(2+m)+68*w+30*w^2+4*w^3)/48;
			end
			if kc <= dim,
				rowi(i) = kr; coli(i) = kc;
				valsA(i) = -2*(1+l)*(-3+14*Z+8*l*Z+n*(-2+4*Z)+m*(-4+8*Z));
				valsB(i) = -2*(1+l)*(14+4*m^2+15*n+n^2+12*l*(1+m+n)+2*m*(9+2*n));
				i = i+1;
			end

			% term no. 30
			if 1+l<=m,
				kc = (285+6*(-1)^(l+m)-3*(-1)^w+12*l^2+48*m+12*m^2+24*l*(4+m)+140*w+42*w^2+4*w^3)/48;
			else
				kc = (237+6*(-1)^(l+m)-3*(-1)^w+12*l^2+96*m+12*m^2+24*l*(2+m)+140*w+42*w^2+4*w^3)/48;
			end
			if kc <= dim,
				rowi(i) = kr; coli(i) = kc;
				valsA(i) = 2*(1+l)*(1+n)*(-1+2*Z);
				valsB(i) = 2*(1+l)*(1+n)*(2+2*m+n);
				i = i+1;
			end

			% term no. 31
			if l<=m,
				kc = (345-6*(-1)^(l+m)-3*(-1)^w+12*l^2+72*m+12*m^2+24*l*(5+m)+140*w+42*w^2+4*w^3)/48;
			else
				kc = (345-6*(-1)^(l+m)-3*(-1)^w+12*l^2+120*m+12*m^2+24*l*(3+m)+140*w+42*w^2+4*w^3)/48;
			end
			if kc <= dim,
				rowi(i) = kr; coli(i) = kc;
				valsA(i) = 4*(1+l)*(1+m)*(-1+2*Z);
				valsB(i) = 4*(1+l)*(1+m)*(2+l+m);
				i = i+1;
			end

			% term no. 32
			if n>=1,
				if 2+l<=m,
					kc = (291-6*(-1)^(l+m)+3*(-1)^w+12*l^2+72*m+12*m^2+24*l*(5+m)+68*w+30*w^2+4*w^3)/48;
				else
					kc = (195-6*(-1)^(l+m)+3*(-1)^w+12*l^2+120*m+12*m^2+24*l*(3+m)+68*w+30*w^2+4*w^3)/48;
				end
				if kc <= dim,
					rowi(i) = kr; coli(i) = kc;
					valsA(i) = 0;
					valsB(i) = 4*(2+3*l+l^2)*n;
					i = i+1;
				end
			end

			% term no. 33
			if 2+l<=m,
				kc = (393-6*(-1)^(l+m)-3*(-1)^w+12*l^2+72*m+12*m^2+24*l*(5+m)+140*w+42*w^2+4*w^3)/48;
			else
				kc = (297-6*(-1)^(l+m)-3*(-1)^w+12*l^2+120*m+12*m^2+24*l*(3+m)+140*w+42*w^2+4*w^3)/48;
			end
			if kc <= dim,
				rowi(i) = kr; coli(i) = kc;
				valsA(i) = 4*(2+3*l+l^2)*Z;
				valsB(i) = 4*(2+3*l+l^2)*(1+m+n);
				i = i+1;
			end


        end
    end
end

matA = sparse(rowi, coli, valsA, dim, dim);
matB = sparse(rowi, coli, valsB, dim, dim);
ev = eigs(matA, matB, 1);

end
