function?CH?=?homogenize（lx?ly?lambda?mu?phi?x）
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%?lx????????=?Unit?cell?length?in?x-direction.
%?ly????????=?Unit?cell?length?in?y-direction.
%?lambda????=?Lame‘s?first?parameter?for?both?materials.?Two?entries.
%?mu????????=?Lame‘s?second?parameter?for?both?materials.?Two?entries.
%?phi???????=?Angle?between?horizontal?and?vertical?cell?wall.?Degrees
%?x?????????=?Material?indicator?matrix.?Size?used?to?determine?nelx/nely
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%?INITIALIZE
%?Deduce?discretization
[nely?nelx]?=?size（x）;
%?Stiffness?matrix?consists?of?two?parts?one?belonging?to?lambda?and
%?one?belonging?to?mu.?Same?goes?for?load?vector
dx?=?lx/nelx;?dy?=?ly/nely;
nel?=?nelx*nely;
[keLambda?keMu?feLambda?feMu]?=?elementMatVec（dx/2?dy/2?phi）;
%?Node?numbers?and?element?degrees?of?freedom?for?full?（not?periodic）?mesh
nodenrs?=?reshape（1:（1+nelx）*（1+nely）1+nely1+nelx）;
edofVec?=?reshape（2*nodenrs（1:end-11:end-1）+1nel1）;
edofMat?=?repmat（edofVec18）+repmat（[0?1?2*nely+[2?3?0?1]?-2?-1]nel1）;
%%?IMPOSE?PERIODIC?BOUNDARY?CONDITIONS
%?Use?original?edofMat?to?index?into?list?with?the?periodic?dofs
nn?=?（nelx+1）*（nely+1）;?%?Total?number?of?nodes
nnP?=?（nelx）*（nely）;????%?Total?number?of?unique?nodes
nnPArray?=?reshape（1:nnP?nely?nelx）;
%?Extend?with?a?mirror?of?the?top?border
nnPArray（end+1:）?=?nnPArray（1:）;
%?Extend?with?a?mirror?of?the?left?border
nnPArray（:end+1）?=?nnPArray（:1）;
%?Make?a?vector?into?which?we?can?index?using?edofMat:
dofVector?=?zeros（2*nn?1）;
dofVector（1:2:end）?=?2*nnPArray（:）-1;
dofVector（2:2:end）?=?2*nnPArray（:）;
edofMat?=?dofVector（edofMat）;
ndof?=?2*nnP;?%?Number?of?dofs
%%?ASSEMBLE?STIFFNESS?MATRIX
%?Indexing?vectors
iK?=?kron（edofMatones（81））‘;
jK?=?kron（edofMatones（18））‘;
%?Material?properties?in?the?different?elements
lambda?=?lambda（1）*（x==1）?+?lambda（2）*（x==2）;
mu?????=?mu（1）*（x==1）?+?mu（2）*（x==2）;
%?The?corresponding?stiffness?matrix?entries
sK?=?keLambda（:）*lambda（:）.‘?+?keMu（:）*mu（:）.‘;
K??=?sparse（iK（:）?jK（:）?sK（:）?ndof?ndof）;?
%%?LOAD?VECTORS?AND?SOLUTION
%?Assembly?three?load?cases?corresponding?to?the?three?strain?cases
sF?=?feLambda（:）*lambda（:）.‘+feMu（:）*mu（:）.‘;
iF?=?repmat（edofMat‘31）;
jF?=?[ones（8nel）;?2*ones（8nel）;?3*ones（8nel）];
F??=?sparse（iF（:）?jF（:）?sF（:）?ndof?3）;
%?Solve?（remember?to?constrain?one?node）
chi（3:ndof:）?=?K（3:ndof3:ndof）\F（3:ndof:）;
%%?HOMOGENIZATION
%?The?displacement?vectors?corresponding?to?