tStart = cputime;
%set side length (lamda), level of approximation (n), 
% and tolerance (tol) if side length is irrational
lamda = 2;
n = 4;
tol = 1e-2;
%gets rational approximation of side length if side length irrational
[N, D] = rat(lamda, tol);
lam = N/D;
indices = subdivide(lam, n);
%form matrix data of holomorphicity + periodicity equations
[row, col, data, rc] = matrixdata(indices, lam, n);
rc=rc+1;
%additional initializer data
row1 = [row, [rc, rc+1, rc+2, rc+2, rc+2,rc+3, rc+3, rc+3]];
col1 = [col, [indices("A1w"), indices("A2w"), indices("B1w"), indices("B1b"), indices("sum1"), indices("B2w"), indices("B2b"), indices("sum2")]];
data1 = [data, [1, 1, 1, 1, -2, 1, 1, -2]];

%two column vectors for two systems of equations
b1 = sparse([rc], [1], [1], rc+4, 1);
b2 = sparse([rc+1], [1], [1], rc+4, 1);

numcols = length(indices);
numrows = rc+4;
%create sparse matrix of equations
M1 = sparse(row1, col1, data1, numrows, numcols);

%solve system
x = M1\b1;
x1 = x(numcols-1);
x2 = x(numcols);
y = M1\b2;
y1 = y(numcols-1);
y2 = y(numcols);

tEnd=cputime-tStart

"level = "+n
%output matrix
[x1, x2; y1, y2]
