c言語のpythonへの書き換え

前提・実現したいこと

C言語で以下のように記述したソースコード(#は意図的に消去)をpythonで書き換えたところ実行結果が大きく違うものが出力されてしまいました。pythonでの値が行ごとに同じ値が出力されてしまいます。原因がわかる方お力を貸していただけるとありがたいです。

該当のソースコード

C
1C
2include<stdio.h>
3include<math.h>
4
5define Z_max 300
6define Delta_z 0.005
7
8define Delta_t 2.0e-12
9define N_max 3000
10
11define epsilon0 8.854e-12
12define mu0 1.256e-6
13
14FILE *p_out;
15int i,n;
16
17double Ey[N_max+1][Z_max+1];
18double Hx[N_max+1][Z_max+1];
19double epsilon[Z_max+1],mu[Z_max+1],sigma[Z_max+1];
20
21void shokichi(){
22
23	for(i=0;i<=Z_max;i++){
24		epsilon[i]=epsilon0;
25		mu[i]=mu0;
26		sigma[i]=0.0;
27
28		if((i>20) && (i<80)){
29			Ey[0][i]=sin(2.0*M_PI*(double)i/20.0);
30			Hx[0][i]=-sin(2.0*M_PI*((double)i-0.5)/20.0)/sqrt(mu[1]/epsilon[i]);
31		}
32		else{
33			Ey[0][i]=0.0;
34			Hx[0][i]=0.0;
35		}
36	}
37}
38int main(){
39	n=0;
40	i=0;
41	p_out=fopen("output.dat","w");
42	shokichi();
43	/*FDTD*/
44	for(n=1;n<=N_max;n++){
45		Ey[n][0]=0;
46		Ey[n][Z_max]=0;
47		for(i=1;i<=Z_max-1;i++){
48			Ey[n][i]=(2*epsilon[i]-sigma[i]*Delta_t)/(2*epsilon[i]+sigma[i]*Delta_t)*Ey[n-1][i]+(Delta_t/Delta_z)*(2/(2*epsilon[i]+sigma[i]*Delta_t))*(Hx[n-1][i+1]-Hx[n-1][i]);
49		};
50		for(i=1;i<=Z_max;i++){
51			Hx[n][i]=Hx[n-1][i]+(1/mu[i])*(Delta_t/Delta_z)*(Ey[n][i]-Ey[n][i-1]);
52		};
53	};
54	/*date*/
55	for(n=0;n<=N_max;n++){
56		for(i=0;i<=Z_max;i++){
57				fprintf(p_out,"%e ",Ey[n][i]);
58		}
59		fprintf(p_out,"\n");
60	}
61	fclose(p_out);
62}

Python
1import math
2import numpy as np
3
4z_max:int=300
5DELTA_Z=0.005
6DELTA_T=2.0e-12
7n_max:int=3000
8EPSILON0=8.854e-12
9MU0=1.256e-6
10
11Ey=[[0]*((z_max)+1)]*((n_max)+1)
12Hx=[[0]*(z_max+1)]*(n_max+1)
13epsilon=[0]*(z_max+1)
14mu=[0]*(z_max+1)
15sigma=[0]*(z_max+1)
16
17
18def ftdt():
19    for i in range (0,z_max+1):
20        epsilon[i]=EPSILON0
21        mu[i]=MU0
22        sigma[i]=0.0
23        if(20<i<80):
24            Ey[0][i]=math.sin(2.0*math.pi*i/20.0)
25            Hx[0][i]=-1*(math.sin(2.0*math.pi*((i-0.5)/20.0))/math.sqrt(mu[1]/epsilon[i]))
26        else:
27            Ey[0][i]=0.0
28            Hx[0][i]=0.0
29            
30            
31ftdt()
32
33for n in range (1,n_max+1):
34    Ey[n][0]=0
35    Ey[n][z_max]=0  
36    for j in range (1,z_max):
37        Ey[n][j]=(2*epsilon[j]-sigma[j]*DELTA_T)/(2*epsilon[j]+sigma[j]*DELTA_T)*Ey[n-1][j]+(DELTA_T/DELTA_Z)*(2/(2*epsilon[j]+sigma[j]*DELTA_T))*(Hx[n-1][j+1]-Hx[n-1][j])
38        
39    for k in range (1,z_max+1):
40        Hx[n][k]=Hx[n-1][k]+(1/mu[k])*(DELTA_T/DELTA_Z)*(Ey[n][k]-Ey[n][k-1])
41        
42f=open('out.dat','w')
43
44
45for l in range(0,n_max+1):
46    for s in range(0,z_max+1):
47        f.write(str(Ey[l][s]))
48    f.write('\n')
49f.close()

試したこと

出力ファイルの中身を比較
C言語で出力されたファイルでは予測された通りの数値が出力されましたが、
pythonでのファイルでは、全体で行ごとに同じ値が出力されています。
↓一部抜粋
0-0.0044853492789758180.008157051769084519-0.010363145703701154
0-0.0044853492789758180.008157051769084519-0.010363145703701154

配列Hxを出力したところやはりpythonの方は行ごとに同じ数値が出力されました。
何度もプログラムを見比べ同じ数式を記述したはずなので、原因がわかりません。