Stack smashing detected C++. Solving Laplace equation

I am solving the Laplace equation using C++. I am solving for P, the pressure field. In my code, you can see below, P is a function of both x and y. I declared it as a 2d array of length nx and ny.

The final matrix result, P, I got, is correct. However, when I tried to print this matrix P using "printf", eventhough the numbers are perfect, I got an error like this:

*** stack smashing detected ***
/bin/bash: line 1: 11738 Aborted

Below is my program:

#include <cmath>
#include <stdio.h>
const int nx = 5; // number of elements in x-direction
const int ny = 5; // number of elements in x-direction
const int niter = 100; //# of iterations 

int main(){
    double dx = 2/double((nx-1)); 
    double dy = 1/double((ny-1));
    int xmax = 2; int xmin = 0;
    int ymax = 1; int ymin = 0;
    double p[nx][ny];
    double pn[nx][ny];
    double x[nx];
    double y[ny];

    //populate p with zeros     
    for (int xi = 1; xi <=nx;xi++){
        for (int yi = 1;yi<=ny;yi++){
            p[xi][yi] = 0;          
        }
    }

    //populate x and y

    //X
    for (int xnum = 1; xnum <=nx; xnum++){
        x[1] = 0;
        x[xnum+1] = x[xnum] + dx;
    }

    //Y
    for (int ynum = 1; ynum <=ny; ynum++){
        y[1] = 0;
        y[ynum+1] = y[ynum] + dy;
    }

    //initial condition

    for (int yrange = 1; yrange<=ny;yrange++){
        p[nx][yrange] = y[yrange];
    }

    //SOLVING FOR P     

    for (int iter = 1; iter<=niter; iter++){
        //copy values
        for (int xiter= 1; xiter<=nx;xiter++){
            for (int yiter = 1; yiter<=ny;yiter++){
                pn[xiter][yiter]=p[xiter][yiter];
            }
        }
        //main loop
        for (int i = 2; i<=nx-1;i++){
            for (int j = 2; j<=ny-1;j++){
                p[i][j] = ((pow(dy,2)*(pn[i+1][j]+pn[i-1][j]))+(pow(dx,2)*(p[i][j+1]+pn[i][j-1])))/(2*(pow(dx,2)+pow(dy,2)));
            }
        }

        for (int xrange = 2; xrange<=nx-1;xrange++){
            p[xrange][1] = p[xrange][2];
            p[xrange][ny] = p[xrange][ny-1];
        }

    }

    //Testing matrix

        for (int x = 1; x<=nx;x++){
            for (int y =1; y<=ny;y++){
            printf("%1.3f\t",p[x][y]);
            }
        printf("\n");
        }

return 0;
}

At first, I thought it's the problem of the number of iteration is too high. In other words, niter is too high, so I decide to lower it. Still, I got the same error. Could you please give me help? The output result of matrix P is correct so I don't think that it is the logic when solving for P. Any help is really appreciated, thank you!

Solution

Your loops go till nx which is undefined behaviour because you are accessing out of bounds. Indexing starts at 0 and goes till i < nx.

This:

for (int xi = 1; xi <= nx; xi++){
    for (int yi = 1; yi <= ny; yi++){
        p[xi][yi] = 0;          
    }
}

Should be

for (int xi = 0; xi < nx; xi++){
    for (int yi = 0; yi < ny; yi++){
        p[xi][yi] = 0;          
    }
}

and the same for all the other loops. Your program may seem to go well and may crash at some place that isnt obviously related to the cause of the UB. When you have undefined behaviour anything can happen (including nothing, so actually you can be happy that you got the crash and have a chance to fix it).

PS: if you really want (not sure why, maybe you are a FORTRAN devotee) you can start your loops at 1. However, I would not recommend to do that, because it makes your loops unnecerassily hard to read and error prone. For example, the error in the following loop can be easily overlooked:

    //copy values (!!! WRONG CODE !!! spot the error to win a FORTRAN fanshirt)
    for (int xiter= 1; xiter < nx+1; xiter++){
        for (int yiter = 1; yiter < ny+1; yiter++){
            pn[xiter-1][yiter-1] = p[xiter][yiter-1];
        }
    }

...btw it is not my intention to bash FORTRAN in any way. In fact I love FORTRAN, but it is a different language (actually the only one I know where you start indices at 1).