% Computes the largest eigenvalue in the ortho case.
% To create a table, type, e.g.:
% >>  for Z=1:10, for om=8:11, tab(Z,om-7) = PekerisOrtho(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 = PekerisOrtho(omega, Z)
dim = floor((omega+1)*(omega+3)*(2*omega+1)/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 l=0:floor((w-1)/2), 
        for m=l+1:w-l,
            n = w-l-m;
            kr = (2*w*(w+2)*(2*w-1) - 3*(1-(-1)^w))/48 + l*(m+n)+m;

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

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

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

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

			% term no. 5
			if l>=1,
				if l<=1+m,
					kc = (3-3*(-1)^w-48*n+48*l*(m+n)-4*w-6*w^2+4*w^3)/48;
					sign = 1;
				else
					kc = (-45-3*(-1)^w+48*l*(1+m)+48*m*(-1+n)-4*w-6*w^2+4*w^3)/48;
					sign = -1;
				end
				if kc <= dim && -1+l ~= m,
					rowi(i) = kr; coli(i) = kc;
					valsA(i) = -2*l*sign*(-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))*sign;
					i = i+1;
				end
			end

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

			% term no. 7
			if l>=1,
				if l<=1+m,
					kc = (-93-3*(-1)^w-48*n+48*l*(2+m+n)+20*w+18*w^2+4*w^3)/48;
					sign = 1;
				else
					kc = (-45-3*(-1)^w+48*l*(1+m)+48*m*(1+n)+20*w+18*w^2+4*w^3)/48;
					sign = -1;
				end
				if kc <= dim && -1+l ~= m,
					rowi(i) = kr; coli(i) = kc;
					valsA(i) = 0;
					valsB(i) = 2*l*(2+3*n+n^2)*sign;
					i = i+1;
				end
			end

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

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

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

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

			% term no. 12
			if m>=1,
				if 1+l<=m,
					kc = (-45-3*(-1)^w+48*m+48*l*(-1+m+n)-4*w-6*w^2+4*w^3)/48;
					sign = 1;
				else
					kc = (3-3*(-1)^w+48*l*m+48*(-1+m)*n-4*w-6*w^2+4*w^3)/48;
					sign = -1;
				end
				if kc <= dim && l ~= -1+m,
					rowi(i) = kr; coli(i) = kc;
					valsA(i) = -2*m*sign*(-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))*sign;
					i = i+1;
				end
			end

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

			% term no. 14
			if m>=1,
				if 1+l<=m,
					kc = (-45-3*(-1)^w+48*m+48*l*(1+m+n)+20*w+18*w^2+4*w^3)/48;
					sign = 1;
				else
					kc = (-93-3*(-1)^w-48*n+48*m*(2+l+n)+20*w+18*w^2+4*w^3)/48;
					sign = -1;
				end
				if kc <= dim && l ~= -1+m,
					rowi(i) = kr; coli(i) = kc;
					valsA(i) = 0;
					valsB(i) = 2*m*(2+3*n+n^2)*sign;
					i = i+1;
				end
			end

			% term no. 15
			if n>=2,
				if l<=m,
					kc = (-3+3*(-1)^w+48*m+48*l*(-2+m+n)+20*w-18*w^2+4*w^3)/48;
					sign = 1;
				else
					kc = (-3+3*(-1)^w+48*l*(1+m)+48*m*(-2+n)+20*w-18*w^2+4*w^3)/48;
					sign = -1;
				end
				if kc <= dim && l ~= m,
					rowi(i) = kr; coli(i) = kc;
					valsA(i) = -((-1+n)*n*sign);
					valsB(i) = 0;
					i = i+1;
				end
			end

			% term no. 16
			if n>=1,
				if l<=m,
					kc = (3-3*(-1)^w+48*m+48*l*(-1+m+n)-4*w-6*w^2+4*w^3)/48;
					sign = 1;
				else
					kc = (3-3*(-1)^w+48*l*(1+m)+48*m*(-1+n)-4*w-6*w^2+4*w^3)/48;
					sign = -1;
				end
				if kc <= dim && l ~= m,
					rowi(i) = kr; coli(i) = kc;
					valsA(i) = 4*n*sign*(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)*sign;
					i = i+1;
				end
			end

			% term no. 17
			if l<=m,
				kc = (-3+3*(-1)^w+48*(1+l)*m+48*l*n-4*w+6*w^2+4*w^3)/48;
				sign = 1;
			else
				kc = (-3+3*(-1)^w+48*l*(1+m)+48*m*n-4*w+6*w^2+4*w^3)/48;
				sign = -1;
			end
			if kc <= dim && l ~= m,
				rowi(i) = kr; coli(i) = kc;
				valsA(i) = 2*sign*(-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)))*sign;
				i = i+1;
			end

			% term no. 18
			if l<=m,
				kc = (3-3*(-1)^w+48*m+48*l*(1+m+n)+20*w+18*w^2+4*w^3)/48;
				sign = 1;
			else
				kc = (3-3*(-1)^w+48*l*(1+m)+48*m*(1+n)+20*w+18*w^2+4*w^3)/48;
				sign = -1;
			end
			if kc <= dim && l ~= m,
				rowi(i) = kr; coli(i) = kc;
				valsA(i) = 4*(1+n)*sign*(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))*sign;
				i = i+1;
			end

			% term no. 19
			if l<=m,
				kc = (45+3*(-1)^w+48*m+48*l*(2+m+n)+68*w+30*w^2+4*w^3)/48;
				sign = 1;
			else
				kc = (45+3*(-1)^w+48*l*(1+m)+48*m*(2+n)+68*w+30*w^2+4*w^3)/48;
				sign = -1;
			end
			if kc <= dim && l ~= m,
				rowi(i) = kr; coli(i) = kc;
				valsA(i) = -((2+3*n+n^2)*sign);
				valsB(i) = 0;
				i = i+1;
			end

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

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

			% term no. 22
			if l<=1+m,
				kc = (51-3*(-1)^w+48*m+48*l*(1+m+n)+20*w+18*w^2+4*w^3)/48;
				sign = 1;
			else
				kc = (3-3*(-1)^w+48*l*(2+m)+48*(1+m)*n+20*w+18*w^2+4*w^3)/48;
				sign = -1;
			end
			if kc <= dim && l ~= 1+m,
				rowi(i) = kr; coli(i) = kc;
				valsA(i) = -2*(1+m)*sign*(-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))*sign;
				i = i+1;
			end

			% term no. 23
			if l<=1+m,
				kc = (93+3*(-1)^w+48*m+48*l*(2+m+n)+68*w+30*w^2+4*w^3)/48;
				sign = 1;
			else
				kc = (93+3*(-1)^w+48*l*(2+m)+48*n+48*m*(1+n)+68*w+30*w^2+4*w^3)/48;
				sign = -1;
			end
			if kc <= dim && l ~= 1+m,
				rowi(i) = kr; coli(i) = kc;
				valsA(i) = 2*(1+m)*(1+n)*sign*(-1+2*Z);
				valsB(i) = 2*(1+m)*(1+n)*(2+2*l+n)*sign;
				i = i+1;
			end

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

			% term no. 25
			if l<=2+m,
				kc = (141+3*(-1)^w+48*m+48*l*(2+m+n)+68*w+30*w^2+4*w^3)/48;
				sign = 1;
			else
				kc = (45+3*(-1)^w+48*l*(3+m)+48*(2+m)*n+68*w+30*w^2+4*w^3)/48;
				sign = -1;
			end
			if kc <= dim && l ~= 2+m,
				rowi(i) = kr; coli(i) = kc;
				valsA(i) = 4*(2+3*m+m^2)*sign*Z;
				valsB(i) = 4*(2+3*m+m^2)*(1+l+n)*sign;
				i = i+1;
			end

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

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

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

			% term no. 29
			if 1+l<=m,
				kc = (3-3*(-1)^w+48*(2+l)*m+48*(1+l)*n+20*w+18*w^2+4*w^3)/48;
				sign = 1;
			else
				kc = (51-3*(-1)^w+48*l*(1+m)+48*m*(1+n)+20*w+18*w^2+4*w^3)/48;
				sign = -1;
			end
			if kc <= dim && 1+l ~= m,
				rowi(i) = kr; coli(i) = kc;
				valsA(i) = -2*(1+l)*sign*(-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))*sign;
				i = i+1;
			end

			% term no. 30
			if 1+l<=m,
				kc = (93+3*(-1)^w+96*m+48*n+48*l*(1+m+n)+68*w+30*w^2+4*w^3)/48;
				sign = 1;
			else
				kc = (93+3*(-1)^w+48*l*(1+m)+48*m*(2+n)+68*w+30*w^2+4*w^3)/48;
				sign = -1;
			end
			if kc <= dim && 1+l ~= m,
				rowi(i) = kr; coli(i) = kc;
				valsA(i) = 2*(1+l)*(1+n)*sign*(-1+2*Z);
				valsB(i) = 2*(1+l)*(1+n)*(2+2*m+n)*sign;
				i = i+1;
			end

			% term no. 31
			if l<=m,
				kc = (141+3*(-1)^w+96*m+48*n+48*l*(1+m+n)+68*w+30*w^2+4*w^3)/48;
				sign = 1;
			else
				kc = (141+3*(-1)^w+48*l*(2+m)+48*n+48*m*(1+n)+68*w+30*w^2+4*w^3)/48;
				sign = -1;
			end
			if kc <= dim && 1+l ~= 1+m,
				rowi(i) = kr; coli(i) = kc;
				valsA(i) = 4*(1+l)*(1+m)*sign*(-1+2*Z);
				valsB(i) = 4*(1+l)*(1+m)*(2+l+m)*sign;
				i = i+1;
			end

			% term no. 32
			if n>=1,
				if 2+l<=m,
					kc = (-93-3*(-1)^w+144*m+96*n+48*l*(-1+m+n)+20*w+18*w^2+4*w^3)/48;
					sign = 1;
				else
					kc = (99-3*(-1)^w+48*l*(1+m)+48*m*(1+n)+20*w+18*w^2+4*w^3)/48;
					sign = -1;
				end
				if kc <= dim && 2+l ~= m,
					rowi(i) = kr; coli(i) = kc;
					valsA(i) = 0;
					valsB(i) = 4*(2+3*l+l^2)*n*sign;
					i = i+1;
				end
			end

			% term no. 33
            if 2+l<=m,
				kc = (45+3*(-1)^w+48*(3+l)*m+48*(2+l)*n+68*w+30*w^2+4*w^3)/48;
				sign = 1;
			else
				kc = (141+3*(-1)^w+48*l*(1+m)+48*m*(2+n)+68*w+30*w^2+4*w^3)/48;
				sign = -1;
            end
            if kc <= dim && 2+l ~= m,
				rowi(i) = kr; coli(i) = kc;
				valsA(i) = 4*(2+3*l+l^2)*sign*Z;
				valsB(i) = 4*(2+3*l+l^2)*(1+m+n)*sign;
				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
