I'm very new to c++ programming, and my first Homework is to implement some iterations for linear system.. is it too much to ask if you could check what's wrong with my code?

here's the code:

my target is to solve a predefined 3500x3500 matrix using jacobi iteration... the code works perfectly for a simple 3x3 matrix (correct answer up to 5 decimal places) however, problem arises when the matrix gets larger (up to 120x120, the matrix still outputs an answer)... iteration always end up with "-1.#IND" error...please please...

any comment or help would be greatly appreciated..many thanks

#include <iostream.h>
#include <math.h>
#define n 121 //number of equations (ok up to 120) [target n=3500]

float a[n][n];
int main()
{
int i, j;
float b[n], h = 1.0/(float)(n+1);
//left hand side matrix of the equation ax=b
for (i=1;i<=n-1;i++)
for (j=1;j<=n-1;j++)
a[i][j] = 0.0;
a[1][1] = -(2.0 + 4.0*h*h);
a[1][2] = 1.0 ;
for (i=2;i<=n-1;i++)
{
a[i][i-1] = 1.0 ;
a[i][i] = -(2.0 + 4.0*h*h);
a[i][i+1] = 1.0 ;
}
a[n][n-1] = 1.0;
a[n][n] = -(2.0 + 4.0*h*h);

//right hand side of the vector of the equation ax=b
for(i=1;i<=n-1;i++)
b[i] = h*h;
b[n] = h*h - 1.0;

//calculation using Jacobi Iteration
double pold[n], pnew[n], tol = 0.000001; //tol = set as tolerance
//pold[]=vector of old values
//pnew[]=vector of new results
int counter = 0, s;
s = 0;
//initialize
for (i=1;i<=n;i++)
pnew[i] = pold[i] = 0.0;
//computation loop
while (s!=1)
{
for(i=1;i<=n;i++)
{
double sum = b[i];
for(j=1;j<=i-1;j++)
sum -= a[i][j]*pold[j];
for(j=i+1;j<=n;j++)
sum -= a[i][j]*pold[j];
pnew[i] = sum/(a[i][i]);
if (fabs(pnew[i]-pold[i])<=tol)
s = 1;
}
counter++;
if (counter%1000 == 0) //indication if program still running
cout << "|";
for (i=1;i<=n;i++)
pold[i] = pnew[i];
}
//print output
cout<<"\n";
for (i=1;i<=n;i++)
cout<<pnew[i]<<endl;
cout<<"\nno. of iterations = "<<counter<<endl;
cin>>i; //system 'PAUSE'
return 0;
}