Thank you DaWei. I apologize for all the inconvenience and unpleasant feelings, but there is absolutely no intention of being rude...

Here's code:

Code:

#include <stdio.h>
#include <math.h>
//double pow(double x, double y);
//#include <grid.h>
//#include <stddev.h>
//#include "ranlib.h"
#define M 100
#define N 200
#define P 100
#define Q 400000
#define X 200 //width
#define Y 400 // height
#define T 500
#define K 0.6
#define MU 0.00089 //viscosity water
#define PHm 0.68
#define n0 20
#define ph0 0.1
#define C 2147483647 // 2^31-1
#define A 16807 //7^5
main()
{
int i, j, k, l, re, rez, rep, p2;
double bin, pi, delx, dely, delt, s, vx1, vx2, vy1, vy2, a, vv, theta, zeta, psi, nn, gau;
double t[P], xx[M+1], yy[N+1], no[M][N][P], D[P][Q], v[M+1][N+1], gammadot[M][N], phi[M][N][P], yita[M][N][P];
double bex[Q][P], bey[Q][P];
//M separations in X direction
//N separations in Y direction
//P no. time intervals
//Q no. particles
pi=4.0*atan(1.0);
re=1; // controls initial theta
rez=16; // controls initial positions of particles
rep=7; // same as above
p2=33; // controls Gaussian distribution
theta=2*pi*re/C;
a=.1;
vv=3*pi*a*a*a/4;
delx=X/M;
dely=Y/N;
delt=T/P;
bin=0;
for(i=0;i<=P-1;i++)
{
t[i]=delt*i;
}
for(i=0;i<=M;i++)
{
xx[i]=i*delx;
for(j=0;j<=N;j++)
{
yy[j]=j*dely;
v[i][j]=0;
for(k=1;k<=5;k++) //the first 5 terms of FFT
{
s=2*k-1;
v[i][j]=v[i][j]+1/s/s/s*(1-cosh(s*pi*(xx[i]-X/2)/Y)/cosh(s*pi*X/(2*Y)))*cos(s*pi*(yy[j]-Y/2)*X/(2*Y*Y));
}
}
}

The debugger targets the 8th last line: v[i][j]=0 that caused the segmentation fault. Also the 5th last line: s=2*k-1 will cause the same error if I disable v[i][j]=0

Thank you for your time and hopefully someone can help out.