## 3 Matlab Implementation
### Iterative Solution (Lines 54-137)
## Abbreviations

The following abbreviations are used in this manuscript:
```
1%%LOOP PARAMETERS
2 clear;
3 NumIter=300;delta=1e-9;
4 nelx=200;nely=100;
5 penal=3;rmin=5;NMaterial=3;
6 NCONN=NMaterial;NDV=nelx*nely*NMaterial;NDV3=nelx*nely;
7 movelimit=0.1;beta=1;
8%%MATERIAL PROPERTIES
9 par.E(1)=1;par.E(2)=2;par.E(3)=5;
10 nu=0.3;
11%%PREPAREFINITEELEMENT ANALYSIS
12A11=[123-6-3;31230;-6312-3;-30-312];
13A12=[-6-303;-3-6-3-6;0-3-63;3-63-6];
14B11=[-43-29;3-4-94;-2-9-4-3;94-3-4];
15B12=[2-34-9;-329-2;4923;-9-232];
16 lk=1/(1-nu^2)/24*([A11A12;A12^A11])=nu*[B11B12;B12^B11]);
17 moderns=reshape(1:(1+nelx)*(1+nely),1+nely);
18 edofVec=reshape(2*nodemrs(1:end-1,1:end-1)+1,nelx*nely,1);
19 edofMat=repmat(edofVec,1,8)+repmat([012*nely+[2301]-2-1],nelx*nely,1);
20 iK=reshape(kron(edofMat,ones(8,1))',64*nelx*nely,1);
21 jK=reshape(kron(edofMat,ones(1,8))',64*nelx*nely,1);
22%%PREPAREFILTER
23 iH=ones(nelx*nely*(2*(ceil(rmin)-1)+1)^2,1);
24 jH=ones(size(iH));
25 sH=zeros(size(iH));
26k=0;
27 fori1=1:nelx
28 forj1=1:nely
29e1=(i1-1)*nely+j1;
30 fori2=max(i1-(ceil(rmin)-1),1):min(i1+(ceil(rmin)-1),nelx)
31 forj2=max(j1-(ceil(rmin)-1),1):min(j1+(ceil(rmin)-1),nely)
32e2=(i2-1)*nely+j2;
33k=k+1;
34 iH(k)=e1;
35 jH(k)=e2;
36 sH(k)=max(0,rmin-sqrt((i1-i2)^2+(j1-j2)^2));
37 end
38 end
39 end
40 end
41H=sparse(iH,jH,sH);
42 Hs=sum(H,2);
43%%INITIALIZE ITERATION
44 loopbeta=0;45 x=ones(NDV,1)*0.5;
46 for i=0:NMaterial-1
47 xTilde(1:nely,1:nelx,i+1)=reshape(x(1+i*NDV3:NDV3+i*NDV3),nely,nelx);
48 end
49 xPhys = (tanh(beta*0.5) + tanh(beta*(xTilde-0.5)))/(tanh(beta*0.5) + tanh(beta*(1-0.5)));
50 % MMA PARAMETTER INITIALIZE
51 a0=0; a=zeros(NCON,1);cc=1.0e6*ones(NCON,1); d=0*ones(NCON,1);
52 xmin=1e-3*ones(NDV,1); xmax=1*ones(NDV,1);
53 xold=x; xolder=xold;low=0; unp=1;
54 %% START ITERATION
55 for iter=1:NumIter
56 loopbeta=loopbeta+1;
57 %% FE-ANALYSIS
58 [KE0] = lk; p=6;
59 for i=1:NMaterial
60 KE(:,:,i)=par.E(i)*KE0;
61 end
62 KE_test=zeros(64,nely*nelx); g=zeros(NMaterial,1);
63 r1=xPhys(:,:,1); r2=xPhys(:,:,2); r3=xPhys(:,:,3);
64 for j=1:NMaterial
65 post_x(1:nely,1:nelx,j)=(xPhys(:,:,j). *(r1.^p + r2.^p + r3.^p ).^(1/p))./(r1 + r2 + r3 + delta);
66 dr(1:nely,1:nelx,j)=(post_x(1:nely,1:nelx,j)).^penal;
67 KE_test=reshape(KE(:,:,j),64,1)*reshape(dr(1:nely,1:nelx,j),1,nely*nelx)+KE_test;
68 g(j)=sum((sum(xPhys(:,:,j)))^));
69 end
70 sK = reshape(KE_test,64*nelx*nely,1);
71 K = sparse(iK,jK,sK); K = (K*K')/2;
72 %% DEFINE LDADS AND SUPPORTS (HALF MBB-BEAM)
73 F = sparse(2*(nely+1)*(nelx+1),1); U = zeros(2*(nely+1)*(nelx+1),1);
74 F(2,1) = -0.1;
75 fixeddofs = union([1:2:2*(nely+1)],[2*(nelx+1)*(nely+1)]);
76 alldofs = [1:2*(nely+1)*(nelx+1)];
77 freedofs = setdiff(alldofs,fixeddofs);
78 %% SOLVING
79 U(freedofs,:) = K(freedofs,freedofs)\F(freedofs,:);
80 U(fixeddofs,:)= 0;
81 c=F'*U;
82 %% SENSITIVITY ANALYSIS
83 dc=zeros(nely,nelx,NMaterial);
84 dgdx_test=zeros(NMaterial,nely*nelx*NMaterial);
85 KES=reshape(sum(U(edofMat)*KE0. *U(edofMat),2),nely,nelx);
86 for n=1:NMaterial
87 for m=1:NMaterial
88 if m==n
89 ddmrn(1:nely,1:nelx,m)=penal. *post_x(1:nelx,m).^(penal -1).*... 90 ((r1.^p + r2.^p + r3.^p).^(1/p)./(r1 + r2 + r3 + delta) -(xPhys(:,:,m). * (r1.^p + r2.^p + r3.^p).^(1/p))./(r1 + r2 + r3 + delta).^2 + (xPhys(:,:,m). *xPhys(:,:,n).^(p -1). *(r1.^p + r2.^p + r3.^p). ^(1/p -1))./(r1 + r2 + r3 + delta));
91 mass_dmdmrn=1;
92 else
93 ddmrn(1:nely,1:nelx,m)=penal. *post_x(1:nely,1:nelx,m).^(penal -1).*... ((xPhys(:,:,m). *xPhys(:,:,n).^(p -1). *(r1.^p + r2.^p + r3.^p). ^(1/p -1))./(r1 + r2 + r3 + delta)- (xPhys(:,:,m). *(r1.^p + r2.^p + r3.^p). \"(1/p))./(r1 + r2 + r3 + delta).^2));
95 mass_dmdrn=0;
96 end
97 dgdx_test(m,1+(n-1)*NDV3:NDV3+NDV3*(n-1))=repmat(mass_dmdrn,1,nelx*nely);
98 end
99 dc(:,:,n)=-dmdrn(:,:,1).*KES.*par.E(1)-dmdrn(:,:,2).*KES.*par.E(2)-dmdrn(:,:,3).*KES. *par.E(3);
100 end
101 %% FILTERING/MODIFICATION OF SENSITIVITIES
102 filtereddc=[];dfdx=[];
103 for i=1:NMaterial
104 dxx = beta* (1-tanh(beta*(xTilde-0.5)).