 function M = matrix2(g, n,Nlam, Dlam, Nmu, Dmu )
%MATRIX gives the sparse matrix encoding
%the equations coming from the Jenkins-Strebel representative of genus g
%and approximation level n


%%HOLOMORPHICITY
%%build entries of the sparse matrix for the holomorphicity conditions
eqstart = 1;
numeqs = (3 + (Nmu/Dmu))*(Dlam*Dmu)*3^n * (Nlam * Dmu)*3^n; %number of equations for holomorphicity conditions
eqend = eqstart + numeqs - 1;
%rows: build four 1s for each of four variables involved
hrow = repmat(eqstart:eqend, [4 1]);
hrow = hrow(:)';
%columns: variables matrix to match the holomorphicity equations
hcol = zeros(1, numeqs*4);
count = 1;
for i = 0:3^n*(3 + (Nmu/Dmu))*(Dlam*Dmu)-1 % Equation: i (x_{i+1, j+1} - x_{i,j}) - (x_{i,j+1} -x_{i+1,j}) = 0 
    for j = 0: (Nlam * Dmu)*3^n-1
        hcol(count) = varindex2(g, n, "NA", i+1, j+1,Nlam, Dlam, Nmu, Dmu );
        hcol(count+1) = varindex2(g, n, "NA", i, j,Nlam, Dlam, Nmu, Dmu );
        hcol(count+2) = varindex2(g, n, "NA", i, j+1,Nlam, Dlam, Nmu, Dmu );
        hcol(count+3) = varindex2(g, n, "NA", i+1, j,Nlam, Dlam, Nmu, Dmu );
        count = count +4;
    end
end
%data: of coefficients of holomorphic equation
hdat = repmat([1i, -1i, -1, 1], 1, numeqs);

%%A1 PERIODS
%build sparse matrix entries corresponding to A periods less than g
eqstart = eqend+1;
numeqs = ((Dlam*Dmu)*3^n+1); %number of equations for first A1 periods
eqend = eqstart + numeqs-1;
%row
aperrow = repmat(eqstart:eqend, [3 1]); %create row vector which repeats eq number
aperrow = aperrow(:)';
%columns
apercol = zeros(1, numeqs*3);
count = 1;

for i=0:(Dlam*Dmu)*3^n %A_{1} = top - bottom
    apercol(count) = varindex2(g, n, "A", 1, mod(i, 2),Nlam, Dlam, Nmu, Dmu );
    apercol(count+1) = varindex2(g, n, "NA", i+(Nmu * Dlam)*3^n, (Nlam*Dmu)*3^n,Nlam, Dlam, Nmu, Dmu );
    apercol(count+2) = varindex2(g, n, "NA", i, 0 ,Nlam, Dlam, Nmu, Dmu );
    count = count +3;
end
%data
aperdat = repmat([-1, 1, -1], 1, numeqs);


%%A2 PERIODS (0 mod 4 SIDES)
eqstart = eqend+1;
numeqs =((Dlam*Dmu)*3^n+1); %sides identified by A_2 period side length 1
eqend = eqstart + numeqs-1;
%rows
agperrow = repmat(eqstart:eqend, [3 1]); %create row vector which repeats eq number
agperrow = agperrow(:)';
%coloumns
agpercol = zeros(1, numeqs*3);
count = 1;


for i = (Dlam*Dmu + Dlam*Nmu)*3^n: (2*Dlam*Dmu + Dlam*Nmu)*3^n  %%%A2 = bottom - top
    agpercol(count) = varindex2(g, n, "A", 2, mod(i, 2),Nlam, Dlam, Nmu, Dmu );
    agpercol(count+1) = varindex2(g, n, "NA", i, 0,Nlam, Dlam, Nmu, Dmu ); 
    agpercol(count+2) = varindex2(g, n, "NA", i+(Dlam*Dmu)*3^n, (Nlam*Dmu)*3^n,Nlam, Dlam, Nmu, Dmu );
    count = count+3;
end

%data
agperdat = repmat([-1, 1, -1], 1, numeqs);


%%B1 PERIOD
eqstart = eqend+1;
numeqs = (Dlam * Nmu)*3^n+1; %length of mu side B periods
eqend = eqstart + numeqs-1;
%row
bperrow = repmat(eqstart:eqend, [3 1]); %create row vector which repeats eq number
bperrow = bperrow(:)';
%columns
bpercol = zeros(1, numeqs*3);
count = 1;

for i=0 :(Dlam*Nmu)*3^n %B_{1} = top - bottom
    bpercol(count) = varindex2(g, n, "B", 1, mod(i+1, 2),Nlam, Dlam, Nmu, Dmu );
    bpercol(count+1) = varindex2(g, n, "NA", i, (Nlam*Dmu)*3^n,Nlam, Dlam, Nmu, Dmu ); 
    bpercol(count+2) = varindex2(g, n, "NA", i+(Dlam*Dmu)*3^n, 0,Nlam, Dlam, Nmu, Dmu );
    count = count+3;
end
%data
bperdat = repmat([-1, 1, -1], 1, numeqs);

%%Periodicity for 3 side %%%B2 = bottom - top
eqstart = eqend+1;
numeqs = (Dlam * Dmu)* 3^n + 1;
eqend = eqstart + numeqs-1;
%row
bgperrow = repmat(eqstart:eqend, [3 1]); %create row vector which repeats eq number
bgperrow = bgperrow(:)';
%columns
bgpercol = zeros(1, numeqs*3);
count = 1;
for i=((Dlam * Dmu) + (Dlam * Nmu))*3^n : (2*(Dlam * Dmu) + (Dlam * Nmu))*3^n%B_{2} = bottom - top
    bgpercol(count) = varindex2(g, n, "B",2, mod(i+1, 2),Nlam, Dlam, Nmu, Dmu );
    bgpercol(count+1) = varindex2(g, n, "NA", i+(Dlam*Dmu)*3^n, 0,Nlam, Dlam, Nmu, Dmu );
    bgpercol(count+2) = varindex2(g, n, "NA", i, (Nlam*Dmu)*3^n,Nlam, Dlam, Nmu, Dmu );
    count = count+3;
end
%data
bgperdat = repmat([-1, 1, -1], 1, numeqs);


%%ZERO SIDES
eqstart = eqend+1;
numeqs = (Nlam * Dmu)*3^n+1;
eqend = eqstart + numeqs-1;
%row
zrow = repmat(eqstart:eqend, [(2+2*g) 1]); %create row vector which repeats eq number
zrow = zrow(:)';
%columns
zcol = zeros(1, numeqs*(2+2*g));
count = 1;
for j = 0:(Nlam*Dmu)*3^n %Oright- 0left = -B1 + A1 + B2 - A2) 
    zcol(count) = varindex2(g, n, "NA", ((3*Dlam*Dmu )+(Dlam*Nmu))*3^n, j,Nlam, Dlam, Nmu, Dmu ); 
    zcol(count+1) = varindex2(g, n, "NA", 0, j,Nlam, Dlam, Nmu, Dmu );
    zcol(count+2) = varindex2(g, n, "B", 1, mod(j, 2),Nlam, Dlam, Nmu, Dmu );
    zcol(count+3) = varindex2(g, n, "A", 1, mod(j, 2),Nlam, Dlam, Nmu, Dmu );
    zcol(count+4) = varindex2(g, n, "B", 2, mod(j, 2),Nlam, Dlam, Nmu, Dmu );
    zcol(count+5) = varindex2(g, n, "A", g, mod(j, 2),Nlam, Dlam, Nmu, Dmu );
    count = count+2+(2*g);
end
%data: create coefficients 
repeater = ones(1, 6);
repeater(1,1) =1; %0right
repeater(1,2) =-1; %0left
repeater(1,3) = 1; %B1
repeater(1,4) =-1; %-A1
repeater(1,5) = -1; %-B2
repeater(1,6) = 1; %A2
zdat = repmat(repeater, 1, numeqs); 

%%SUMS 
%variable sum is given by average of black and white B values
eqstart = eqend+1;
numeqs = 2;
eqend = eqstart + numeqs -1;
%row
srow = repmat(eqstart:eqend, [3 1]); %create row vector which repeats eq number
srow = srow(:)';
%columns
scol = zeros(1, numeqs*3);
count = 1;
for i = 1:g
    scol(count) = varindex2(g, n, "B", i, 0,Nlam, Dlam, Nmu, Dmu );
    scol(count+1) = varindex2(g, n, "B", i, 1,Nlam, Dlam, Nmu, Dmu );
    scol(count+2) = varindex2(g, n, "sum", i, 0,Nlam, Dlam, Nmu, Dmu );
    count = count+3;
end
%data
sdat = repmat([-1, -1, 2], 1, numeqs);


%%INITIALIZERS
%Ab's=Aw's
eqstart = eqend+1;
%rows
erow = repmat(eqstart:eqstart+g-1, [2 1]); %create row vector which repeats eq number
erow = erow(:)';
ecol = zeros(1, g*2);
count = 1;
%columns
for i = 1:g
    ecol(count) = varindex2(g, n, "A", i, 0,Nlam, Dlam, Nmu, Dmu );
    ecol(count+1) = varindex2(g, n, "A", i, 1,Nlam, Dlam, Nmu, Dmu );
    count = count+2;
end
%data
edat = repmat([1, -1], 1, g);

%Canonical basis so A_{ij} = delta_{ij}
eqend = eqstart+g-1;
%rows
rowinit = eqend+1:eqend+g+2;
%columns
colinit = zeros(1, g);
for i = 1:g
    colinit(i) = varindex2(g, n, "A", i, 0,Nlam, Dlam, Nmu, Dmu );
end
colinit = [varindex2(g, n, "NA", 0, 0,Nlam, Dlam, Nmu, Dmu ), varindex2(g, n, "NA", 1, 0,Nlam, Dlam, Nmu, Dmu ), colinit];
%data
datainit = ones(1, g+2);


%make final row, columns, data to construct the sparse matrix
row = [hrow, aperrow, agperrow, bperrow, bgperrow, zrow, srow, erow, rowinit];
col = [hcol, apercol, agpercol,  bpercol, bgpercol, zcol, scol, ecol, colinit];
data = [hdat, aperdat, agperdat,  bperdat, bgperdat, zdat, sdat, edat, datainit];
M = sparse(row, col, data);
end