連立一次方程式を用いた二次元の電場と磁場の波のシミュレーションのプログラムについてお願いします。

[ヴぁヴぁ]

連立一次方程式を用いた二次元の電場と磁場の波のシミュレーションのプログラムについてなのですが、
連立一次方程式の計算(Gauss消去法)をサブルーチンで計算させるプログラムを作りました。
その結果コンパイルはでき走らせるまでは行くのですが数値が-1.#QOと出てしまい求める事ができません。
どなたかわかる方いらっしゃいましたらアドバイス等よろしくお願いします。
(hx、hy、dz、ezの計算は当方で考えた計算ですので、プログラムとしてのミスではなく計算のミスと思われる場合はそのように言っていただけると幸いです。)

nstepsの数字は任意ですので何でもかまいません。(当方は5で走らせていました。)

以下がプログラムになります。

#include <math.h>
#include <stdlib.h>
#include <stdio.h>

#define IE 60
#define JE 60

void lude(float at[][3], float b[], int n )
{
int k;
for(k=1;k<=n;k++)
{
at[k][1]=1.;
at[k][2]=-2.25;
at[k][3]=1.;
}

for(k=1;k<n;k++)
{

b[k]=b[k]/at[k][2];
at[k][3]=at[k][3]/at[k][2];
b[k+1]=b[k+1]-at[k+1][1]*b[k];
at[k+1][2]=at[k+1][2]-at[k+1][1]*at[k][3];
}

b[n]=b[n]/at[n][2];

for(k=n;k>=1;k--)
{
b[k]=b[k]-at[k][3]*b[k+1];
}
}

main()
{
float ga[IE][JE],dz[IE][JE],ez[IE][JE],hx[IE][JE],hy[IE][JE];
float at[JE][3],b[JE];
int l,n,i,j,ic,jc,nsteps,k;
float ddx,dt,T,epsz,pi,epsilon,sigma,eaf;
float t0,spread,pulse;
FILE *fp, *fopen();

for(j=0;j<JE;j++)
{
ic=IE/2;
jc=JE/2;
ddx=.01; /* Cell size */
dt=ddx/6e8; /* Time steps */
epsz = 8.85419e-12;
pi=3.141592653;

for(j=0;j<JE;j++)
{
printf("%2d",j);
for(i=0;i<IE;i++)
{
dz[j]=0.;
ez[j]=0.;
hx[j]=0.;
hy[j]=0.;
ga[j]=1.0;
printf("%5.2f ",ga[j]);
}
printf("\n");
}
t0=20.0;
spread=6.0;
T=0.;
nsteps=1.;

printf("nsteps-->");
scanf("%d",&nsteps);
printf("%d\n",nsteps);
n = 0;

for(n=1;n<=nsteps;n++)
{
T=T+1;

/* --- start of the Main FDTD loop --- */


/* Calculate the Hx field */
for(i=0;i<IE-1;i++)
{
for(j=0;j<=JE-1;j++)
{

b[j]=-16.0*hx[j]+4.0*(ez[j+1]-ez[j])+(hy[j+1]-hy[i-1][j]-hy[i][j]+hy[i-1][j]);

}
lude( at, b, JE);

for(j=0;j<JE-1 ;j++)
hx[i][j]=b[j];
}

/* Calculate the Hy field */
for(j=0;j<JE-1;j++)
{
for(i=0;i<=IE-1;i++)
{
hy[i][j]=hy[i][j]+.5*(ez[i+1][j]-ez[i][j]);
}
}

/* Calculate the Dz field */
for(j=1;j<IE;j++)
{
for(i=1;i<IE;i++)
{
dz[i][j]=dz[i][j]+.5*(hy[i][j]-hy[i-1][j]-hx[i][j]+hx[i][j-1]);
}
}
/* Put a Gaussian Pulse the middle */
pulse=exp(-.5*(pow((t0-T)/spread,2.0)));
dz[ic][jc]=pulse;

/* Calculate the Ez field */
for(j=1;j<JE;j++)
{
for(i=1;i<IE;i++)
{
ez[i][j]=ga[i][j]*dz[i][j];
}
}

/* Calculate the Hx field */
for(j=0;j<JE-1;j++)
{
for(i=0;i<=IE-1;i++)
{
hx[i][j]=hx[i][j]-.5*(ez[i][j+1]-ez[i][j]);
}
}

/* Calculate the Hy field */
for(i=1;i<=IE-1;i++)
{
for(j=0;j<JE-1;j++)
{

b[j]=-16*hy[i][j]-4*(ez[i+1][j]-ez[i][j])+(hx[i+1][j]-hx[i+1][j-1]-hx[i][j]+hx[i][j-1]);
}

lude( at, b, JE);

for( j=0; j<JE-1; j++)
hy[i][j]=b[j];
}

/* Calculate the Dz field */
for(j=1;j<IE;j++)
{
for(i=1;i<IE;i++)
{
dz[i][j]=dz[i][j]+.5*(hy[i][j]-hy[i-1][j]-hx[i][j]+hx[i][j-1]);
}
}

/* Put a Gaussian Pulse the middle */
pulse=exp(-.5*(pow((t0-T)/spread,2.0)));
dz[ic][jc]=pulse;

/* Calculate the Ez field */
for(j=1;j<JE;j++)
{
for(i=1;i<IE;i++)
{
ez[i][j]=ez[i][j]+(ga[i][j]*dz[i][j]);
}
}

[ライス]

すみません。
ご質問の内容とは別なのですが、コードを貼り付ける場合はコートタグを使うと便利ですよ。
というより、一応フォーラムルールでもあります。
投稿する際に、
［ｃｏｄｅ］
ここ
［／ｃｏｄｅ］
として、ここと書いているところを削除した後コードを貼リ付けてみてください。
（［ｃｏｄｅ］［／ｃｏｄｅ］は表示されるように全角で書いていますが、半角で書いてください。）

コード:

#include <math.h>
#include <stdlib.h>
#include <stdio.h>

#define IE 60
#define JE 60

void lude(float at[][3], float b[], int n )
{
int k;
for(k=1;k<=n;k++)
{
at[k][1]=1.;
at[k][2]=-2.25;
at[k][3]=1.;
}

for(k=1;k<n;k++)
{

b[k]=b[k]/at[k][2];
at[k][3]=at[k][3]/at[k][2];
b[k+1]=b[k+1]-at[k+1][1]*b[k];
at[k+1][2]=at[k+1][2]-at[k+1][1]*at[k][3];
}

b[n]=b[n]/at[n][2];

for(k=n;k>=1;k--)
{
b[k]=b[k]-at[k][3]*b[k+1];
}
}

main()
{
float ga[IE][JE],dz[IE][JE],ez[IE][JE],hx[IE][JE],hy[IE][JE];
float at[JE][3],b[JE];
int l,n,i,j,ic,jc,nsteps,k;
float ddx,dt,T,epsz,pi,epsilon,sigma,eaf;
float t0,spread,pulse;
FILE *fp, *fopen();

for(j=0;j<JE;j++)
{
ic=IE/2;
jc=JE/2;
ddx=.01; /* Cell size */ 
dt=ddx/6e8; /* Time steps */
epsz = 8.85419e-12;
pi=3.141592653;

for(j=0;j<JE;j++)
{
printf("%2d",j);
for(i=0;i<IE;i++)
{
dz[i][j]=0.;
ez[i][j]=0.;
hx[i][j]=0.;
hy[i][j]=0.;
ga[i][j]=1.0;
printf("%5.2f ",ga[i][j]);
}
printf("\n");
}
t0=20.0;
spread=6.0;
T=0.;
nsteps=1.;

printf("nsteps-->");
scanf("%d",&nsteps);
printf("%d\n",nsteps);
n = 0;

for(n=1;n<=nsteps;n++)
{
T=T+1;

/* --- start of the Main FDTD loop --- */


/* Calculate the Hx field */
for(i=0;i<IE-1;i++)
{
for(j=0;j<=JE-1;j++)
{

b[j]=-16.0*hx[i][j]+4.0*(ez[i][j+1]-ez[i][j])+(hy[i][j+1]-hy[i-1][j]-hy[i][j]+hy[i-1][j]);

}
lude( at, b, JE);

for(j=0;j<JE-1 ;j++)
hx[i][j]=b[j];
}

/* Calculate the Hy field */
for(j=0;j<JE-1;j++)
{
for(i=0;i<=IE-1;i++)
{
hy[i][j]=hy[i][j]+.5*(ez[i+1][j]-ez[i][j]);
}
}

/* Calculate the Dz field */
for(j=1;j<IE;j++)
{
for(i=1;i<IE;i++)
{
dz[i][j]=dz[i][j]+.5*(hy[i][j]-hy[i-1][j]-hx[i][j]+hx[i][j-1]);
}
}
/* Put a Gaussian Pulse the middle */
pulse=exp(-.5*(pow((t0-T)/spread,2.0)));
dz[ic][jc]=pulse;

/* Calculate the Ez field */
for(j=1;j<JE;j++)
{
for(i=1;i<IE;i++)
{
ez[i][j]=ga[i][j]*dz[i][j];
}
}

/* Calculate the Hx field */
for(j=0;j<JE-1;j++)
{
for(i=0;i<=IE-1;i++)
{
hx[i][j]=hx[i][j]-.5*(ez[i][j+1]-ez[i][j]);
}
}

/* Calculate the Hy field */
for(i=1;i<=IE-1;i++)
{
for(j=0;j<JE-1;j++)
{

b[j]=-16*hy[i][j]-4*(ez[i+1][j]-ez[i][j])+(hx[i+1][j]-hx[i+1][j-1]-hx[i][j]+hx[i][j-1]);
}

lude( at, b, JE);

for( j=0; j<JE-1; j++) 
hy[i][j]=b[j];
}

/* Calculate the Dz field */
for(j=1;j<IE;j++)
{
for(i=1;i<IE;i++)
{
dz[i][j]=dz[i][j]+.5*(hy[i][j]-hy[i-1][j]-hx[i][j]+hx[i][j-1]);
}
}

/* Put a Gaussian Pulse the middle */
pulse=exp(-.5*(pow((t0-T)/spread,2.0)));
dz[ic][jc]=pulse;

/* Calculate the Ez field */
for(j=1;j<JE;j++)
{
for(i=1;i<IE;i++)
{
ez[i][j]=ez[i][j]+(ga[i][j]*dz[i][j]);
}
}

このように、見やすくなりますよ。

[box]

全体を見たわけではないので何とも言えませんが、少なくとも

＞void lude(float at[][3], float b[], int n )

at[][3]
という定義もしくは宣言をしている状態で、

＞ at[k][3]=1.;

これは大いにまずいですね。配列の範囲外の領域を参照しています。
配列atの第2次元で有効なのは0～2です。

他にもまずい点があるかもしれませんが、とりあえずここだけ。

# 人に見せるソースコードには、適切なコメントを書いてほしい。
# どの変数がどういう意味を持つのか、さっぱりわからない。