关于错误-1.#IND00

[cyyh0227]

代码如下：
#include<stdio.h>
#include<math.h>
double sin(double x);
#define n 60
#define PI 3.1415926
double a=0.02;
double t=1;
double c=0,d=1; //区间【0,1】
double k=(d-c)/n; //空间方向步长
double h=0.01; //时间方向步长
int m=(int)(t/h); //时间方向迭代次数


int main()
{
int i,j;
int flag=0;
double U[n];
double N[n];
double V[n];
double K1[n]={0},K2[n]={0},K3[n]={0},K4[n]={0};
double L[n][n]={0};
double M[n][n]={0};

for(i=0;i<n;i++)
{
U[i]=sin(i*k*PI);
}

//对矩阵L赋初始值
L[0][0]=-2;
L[0][1]=1;
L[n-1][n-2]=1;
L[n-1][n-1]=-2;
for(i=1;i<n-1;i++)
{
L[i][i-1]=1,L[i][i]=-2,L[i][i+1]=1;
}
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
{
L[i][j]=L[i][j]*a/k/k;
}
}
//L赋值完成


//对矩阵M赋初始值
M[0][0]=U[1];
M[n-1][n-1]=(-1)*U[n-2];
for(i=1;i<n-1;i++)
{
M[i][i]=U[i+1]-U[i-1];
}
for(i=0;i<n;i++)
{
M[i][i]=M[i][i]*(-1)*0.5/k;
}
//M赋值完成



while(flag<m)
{
for(i=0;i<n;i++)
{
V[i]=U[i];
}
//数组V1=U

for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
{
K1[i]=K1[i]+L[i][j]*V[j]+M[i][j]*V[j];
}
}
//计算K1

for(i=0;i<n;i++)
{
V[i]=U[i]+0.5*h*K1[i];
}

for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
{
K2[i]=K2[i]+L[i][j]*V[j]+M[i][j]*V[j];
}
}
//计算K2

for(i=0;i<n;i++)
{
V[i]=U[i]+0.5*h*K2[i];
}

for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
{
K3[i]=K3[i]+L[i][j]*V[j]+M[i][j]*V[j];
}
}
//计算K3

for(i=0;i<n;i++)
{
V[i]=U[i]+h*K3[i];
}

for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
{
K4[i]=K4[i]+L[i][j]*V[j]+M[i][j]*V[j];
}
}
//计算K4

for(i=0;i<n;i++)
{
U[i]=U[i]+(K1[i]+2*K2[i]+2*K3[i]+K4[i])*h/6;
}

for(i=0;i<n;i++)
{
K1[i]=0,K2[i]=0,K3[i]=0,K4[i]=0;
}
//再次将K1,K2，K3，K4初始化


M[0][0]=U[1];
M[n-1][n-1]=(-1)*U[n-2];
for(i=1;i<n-1;i++)
{
M[i][i]=U[i+1]-U[i-1];
}
for(i=0;i<n;i++)
{
M[i][i]=M[i][i]*(-1)*0.5/k;
}
//改变M矩阵的数值

flag++;
}

for(i=0;i<n;i++) printf("U[%d]=%lf\n",i,U[i]);


return 0;
}

这是用龙格库塔法求解方程的程序。但是当n>=60的时候就出现-1.#IND00，而n<60的时候结果无误。