## 3 Matlab Implementation
### Iterative Solution (Lines 54-137)  
These lines are the process of iteratively solving the optimization model, which is the critical part of the implementation. The finite element analysis, the computation of the interpolation function, and the mass constraint for each material are carried out in lines 57-71, whereas the sparse global stiffness matrix is assembled in line 71, as shown in Box 1. The formulation of each material \(\psi_{i}\) is stored in matrix post_x, and the constraints are stored in matrix g.  
```
57%%FE-ANALYSIS
58[XEO]=lk;p=6;
59 fori=1:NMaterial
60 KE(:,:,i)=par.E(i)*KEO;
61 end
62 KE_test=zeros(64,nely*nelx);g=zeros(NMaterial,1);
63r1=xPhys(:,:,1);r2=xPhys(:,:,2);r3=xPhys(:,:,3);
64 forj=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;
68g(j)-sum(sum(xPhys(:,:,j))^()));
69 end
70 sK=reshape(KE_test,64*nelx*nely,1);
71K=sparse(iK,jK,sK);K=(K+K^)/2;The boundary conditions and the load vector of the half MBB beam are defined in lines 72-77. The objective function c is computed in lines 79-81, according to the global stiffness matrix K and the load vector F. The sensitivities dc and dgdx_test of the objective function and the volume constraints with respect to the projected density are computed in lines 82-100, and then the sensitivities of the filtered and projected procedure are added in lines 101-113, which are stored in dfdx and dgdx, as shown in Box 2. ```
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 forn=1:NMaterial
87 form=1:NMaterial
88 ifm==m
89 dmdrn(1:nely,1:nelx,m)=penal.*post_x(1:nely,1:nelx,m).^(pamal-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(:,:,m).^(p-1). * (r1.^p+r2.^p+r3.^p).^(1/p-1))./(r1+r2+r3+delta));
91 mass_dmdmr=1;
92 else
93 dmdm(1:nely,1:nelx,m)=penal.*post_x(1:nely,1:nelx,m).^(pamal-1).*... 94(((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_dmdmr=0;
96 end
97 dgdx_test(m,1+(n-1)*NDV3:NDV3+NDV3*(n-1))= repmat(mass_dmdmr,1,nelx=nely);
98 end
99 dc(:,:,n)=-dmdmr(:,:,1).*KES.*par.E(1)-dmdmr(:,:,2).*KES.*par.E(2)-dmdmr(:,:,3).*KES. *par.E(3);
100 end
101%%FILTERING/MODIFICATION OFSENSITIVITIES
102 filtereddc=[];dfdx=[];dfdx=[];
103 fori=1:NMaterial
104 dxx=beta*(1-tanh(beta*(xTilde-0.5)). *tanh(beta*(xTilde-0.5)))/ (tanh(beta*0.5)+tanh(beta*(1-0.5)));
105[filtereddc(:,:,:)]=H*(reshape(dc(:,:,i). *dxx(:,:,i),NDV3,1)./Hs);
106 dfdx=[dfdx;filtereddc(:,:,:)];
107g(i)=(g(i)/NDV3-0.5/3)*1;
108 dgdx1=(dgdx_test(i,:)/NDV3)*1;
109 forj=0:NMaterial-1
110 dgdx2=reshape(dgdx1(:,:+j*NDV3:NDV3+j*NDV3),nelx);
111 dgdx(i,1+j*NDV3:NDV3+j*NDV3)= H*(reshape(dgdx2(:,:). *dxx(:,:,i),NDV3,1)./Hs);
112 end
113 end
```
The objective function and its sensitivity are scaled to balance the objective and constraints in line 115 if needed, and the movement limits of the MMA algorithm are updated in lines 116-118. The MMA optimization is conducted in lines 119-123, and the intermediate results are displayed in line 125. Parameters are updated in lines 127-136 for the next iteration. The projection function threshold beta can be directly changed to a larger value or be gradually increased from a small value. The influence of the two methods will be discussed later. We start with beta = 1 and double its value after 50 iterations here. ### Image Display (Lines 138-145)

The final design is displayed in lines 139-144, and the optimized structures with multiple materials are shown in different colors. Finally, the result is saved in line 145, as shown in Box 3. **Box 3. ** Image Display code.