Matrix Calculation Different between Python and Matlab

I am transferring a finite element code from Matlab to Python.

A problem occurs at the last step when I try to solve the displacement $U = F/K$.
I have checked that the calculated $F$ and $K$ are same for Matlab and Python. However the resulting $U$ are different.

Can anyone suggest what’s the problem here and how can I solve it?

Below are the codes for Matlab and Python.

Matlab

clear all; clc; close all
nelx=80; 
nely=80; 
beta = 1;
eta = 0.47;
volrate=0.20; 
E0 = 1;
Emin = 0.000001;
nu = 0.3;
a = 0.005;
b = 0.005;
h = 0.01;
q = -0.001;
penal = 3;
[KE1]= ElementStiffness(E0,nu,a,b,h);
TM3 = zeros(nely,nelx);
TM4 = zeros(nely+1,nelx+1);
TM3 = TM3+0.01;
% TM4([98:101 199:202 300:303 401:404]) = 0.01;
TM4([1]) = 0.01;
T2 = zeros((nely+1)*(nelx+1),1); 
num = 0;
for i=1:nely+1
    for j =1:nelx+1
        num = num+1;
        T2(num,1) = TM4(j,i);
    end
end
xTilde = repmat(volrate,[nely nelx]);
xPhys = ((tanh(beta*eta)+tanh(beta*(xTilde-eta)))/(tanh(beta*eta)+tanh(beta*(1-eta))));


K=sparse(3*(nelx+1)*(nely+1),3*(nelx+1)*(nely+1)); 
F=sparse(3*(nely+1)*(nelx+1),1);
U1=sparse(3*(nely+1)*(nelx+1),1); 

Fe=q*a*b*[1 b/3 -a/3 1 b/3 a/3 1 -b/3 a/3 1 -b/3 -a/3]'; 

for ely=1:nely 
    for elx=1:nelx 
        n1 = (nely+1)*(elx-1)+ely;
        n2 = (nely+1)* elx +ely;
        edof = [3*n1+1; 3*n1+2; 3*n1+3;3*n2+1;3*n2+2;3*n2+3;3*n2-2;3*n2-1; 3*n2; 3*n1-2;3*n1-1; 3*n1];
        K(edof,edof)=K(edof,edof)+(Emin+(1-Emin)*xPhys(ely,elx)^penal)*KE1; 
        if TM3(ely,elx)>0 
            F(edof)=F(edof)+Fe;
        end
    end 
end 


dd = find(T2 == 0.01); %T2=0-fixed node
fixeddofs = [(dd-1)*3+1 (dd-1)*3+2 (dd-1)*3+3]; 
alldofs   =1:3*(nely+1)*(nelx+1); 
freedofs  =setdiff(alldofs,fixeddofs); 

U1(freedofs,:)=K(freedofs,freedofs)F(freedofs,:); 
U1(fixeddofs,:)=0;
U1matlab = full(U1);
function [KE1]= ElementStiffness(E,nu,a,b,h)
syms x y 
N1  = (1 - x/a)*(1 - y/b)*(2 - x/a - y/b - (x/a)^2 - (y/b)^2)/8;
Nx1 = b*(1 - x/a)*(1 - y/b)*(1 - (y/b)^2)/8;
Ny1 =  - a*(1 - x/a)*(1 - y/b)*(1 - (x/a)^2)/8;
N2  = (1 + x/a)*(1 - y/b)*(2 + x/a - y/b - (x/a)^2 - (y/b)^2)/8;
Nx2 = b*(1 + x/a)*(1 - y/b)*(1 - (y/b)^2)/8;
Ny2 =   a*(1 + x/a)*(1 - y/b)*(1 - (x/a)^2)/8;
N3  = (1 + x/a)*(1 + y/b)*(2 + x/a + y/b - (x/a)^2 - (y/b)^2)/8;
Nx3 = -b*(1 + x/a)*(1 + y/b)*(1 - (y/b)^2)/8;
Ny3 =   a*(1 + x/a)*(1 + y/b)*(1 - (x/a)^2)/8;
N4  = (1 - x/a)*(1 + y/b)*(2 - x/a + y/b - (x/a)^2 - (y/b)^2)/8;
Nx4 = -b*(1 - x/a)*(1 + y/b)*(1 - (y/b)^2)/8;
Ny4 =  - a*(1 - x/a)*(1 + y/b)*(1 - (x/a)^2)/8;
N=[N1 Nx1 Ny1 N2 Nx2 Ny2 N3 Nx3 Ny3 N4 Nx4 Ny4];
B1=diff(N,x,2);
B2=diff(N,y,2);
B3=diff(N,x,y);
B=-[B1;B2;2*B3];
D = (E/(1-nu*nu))*[1, nu, 0 ; nu, 1, 0 ; 0, 0, (1-nu)/2];
BD = transpose(B)*D*B;
KE1=int(int(BD,x,-a,a),y,-b,b)*h^3/12;
KE1=double(KE1);
end

Python

from __future__ import division 
import numpy as np 
from scipy.sparse.linalg import spsolve
from scipy import sparse

nelx=80
nely=80
VF = 0.2
beta = 1
eta = 0.47
penal = 3 
a = 0.005
b = 0.005
q = -0.001
Emin = 1e-6
TM3 = np.zeros((nely,nelx))
TM3 = TM3+0.01
TM4 = np.zeros((nely+1,nelx+1))
TM4[0,0] = 0.01
T2 = np.zeros(((nely+1)*(nelx+1),1))
num = 0
for i in range(nely+1):
    for j in range(nelx+1):       
        T2[num,0] = TM4[j,i]
        num = num+1
KE1 = [[ 9.67032967e-03,  2.23443223e-05, -2.23443223e-05,
    -4.17582418e-03,  5.12820513e-06, -1.95970696e-05,
    -1.31868132e-03,  7.87545788e-06, -7.87545788e-06,
    -4.17582418e-03,  1.95970696e-05, -5.12820513e-06],
   [ 2.23443223e-05,  1.39194139e-07, -2.74725275e-08,
     5.12820513e-06,  4.39560440e-08,  0.00000000e+00,
    -7.87545788e-06,  3.47985348e-08,  0.00000000e+00,
    -1.95970696e-05,  5.67765568e-08,  0.00000000e+00],
   [-2.23443223e-05, -2.74725275e-08,  1.39194139e-07,
     1.95970696e-05,  0.00000000e+00,  5.67765568e-08,
     7.87545788e-06,  0.00000000e+00,  3.47985348e-08,
    -5.12820513e-06,  0.00000000e+00,  4.39560440e-08],
   [-4.17582418e-03,  5.12820513e-06,  1.95970696e-05,
     9.67032967e-03,  2.23443223e-05,  2.23443223e-05,
    -4.17582418e-03,  1.95970696e-05,  5.12820513e-06,
    -1.31868132e-03,  7.87545788e-06,  7.87545788e-06],
   [ 5.12820513e-06,  4.39560440e-08,  0.00000000e+00,
     2.23443223e-05,  1.39194139e-07,  2.74725275e-08,
    -1.95970696e-05,  5.67765568e-08,  0.00000000e+00,
    -7.87545788e-06,  3.47985348e-08,  0.00000000e+00],
   [-1.95970696e-05,  0.00000000e+00,  5.67765568e-08,
     2.23443223e-05,  2.74725275e-08,  1.39194139e-07,
     5.12820513e-06,  0.00000000e+00,  4.39560440e-08,
    -7.87545788e-06,  0.00000000e+00,  3.47985348e-08],
   [-1.31868132e-03, -7.87545788e-06,  7.87545788e-06,
    -4.17582418e-03, -1.95970696e-05,  5.12820513e-06,
     9.67032967e-03, -2.23443223e-05,  2.23443223e-05,
    -4.17582418e-03, -5.12820513e-06,  1.95970696e-05],
   [ 7.87545788e-06,  3.47985348e-08,  0.00000000e+00,
     1.95970696e-05,  5.67765568e-08,  0.00000000e+00,
    -2.23443223e-05,  1.39194139e-07, -2.74725275e-08,
    -5.12820513e-06,  4.39560440e-08,  0.00000000e+00],
   [-7.87545788e-06,  0.00000000e+00,  3.47985348e-08,
     5.12820513e-06,  0.00000000e+00,  4.39560440e-08,
     2.23443223e-05, -2.74725275e-08,  1.39194139e-07,
    -1.95970696e-05,  0.00000000e+00,  5.67765568e-08],
   [-4.17582418e-03, -1.95970696e-05, -5.12820513e-06,
    -1.31868132e-03, -7.87545788e-06, -7.87545788e-06,
    -4.17582418e-03, -5.12820513e-06, -1.95970696e-05,
     9.67032967e-03, -2.23443223e-05, -2.23443223e-05],
   [ 1.95970696e-05,  5.67765568e-08,  0.00000000e+00,
     7.87545788e-06,  3.47985348e-08,  0.00000000e+00,
    -5.12820513e-06,  4.39560440e-08,  0.00000000e+00,
    -2.23443223e-05,  1.39194139e-07,  2.74725275e-08],
   [-5.12820513e-06,  0.00000000e+00,  4.39560440e-08,
     7.87545788e-06,  0.00000000e+00,  3.47985348e-08,
     1.95970696e-05,  0.00000000e+00,  5.67765568e-08,
    -2.23443223e-05,  2.74725275e-08,  1.39194139e-07]]     
KE1 = np.array(KE1)
xTilde = np.ones(nelx*nely)*VF
xPhys = ((np.tanh(beta*eta)+np.tanh(beta*(xTilde-eta)))/(np.tanh(beta*eta)+np.tanh(beta*(1-eta))))

K = np.zeros((3*(nelx+1)*(nely+1),3*(nelx+1)*(nely+1)))
F =np.zeros((3*(nely+1)*(nelx+1),1))
U1 = np.zeros((3*(nely+1)*(nelx+1),1))
Fe=q*a*b*np.mat([1,b/3,-a/3,1,b/3,a/3,1,-b/3,a/3,1,-b/3,-a/3]).T
xPhys=xPhys.reshape(nelx,nely)

for ely in range(nely):
    for elx in range(nelx):
        n1 = (nely+1)*(elx+1-1)+ely+1
        n2 = (nely+1)*(elx+1) +ely+1
        edof = [3*n1+1, 3*n1+2, 3*n1+3,3*n2+1,3*n2+2,3*n2+3,3*n2-2,3*n2-1, 3*n2, 3*n1-2,3*n1-1, 3*n1]
        edof = [m-1 for m in edof]
        indices = np.ix_(edof, edof)
        K[indices]=K[indices]+(Emin+(1-Emin)*xPhys[ely,elx]**penal)*KE1           
        F[edof,:]=F[edof,:]+Fe   
        
dd = np.where(T2 == 0.01)[0]        
fixeddofs = [(dd)*3+1,(dd)*3+2,(dd)*3+3]
fixeddofs = np.array(fixeddofs)-1  
alldofs = range(3*(nely+1)*(nelx+1))
freedofs =np.setdiff1d(alldofs,fixeddofs).tolist()
K = sparse.csc_matrix(K)
indices = np.ix_(freedofs, freedofs)
U1[freedofs,0]=spsolve(K[indices],F[freedofs,0])
U1[fixeddofs,:]=0