Erro: Leapfrog

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


void atualizar_onda(double u_new[],double u_now[],double u_old[],int jmax,double k);

void main()
{

    FILE *arq;
    arq = fopen("erro_leapfrog.dat", "w");

    int i, j, jmax;
    
    //tamanho da corda: jmax-1
    jmax = 50;
    
    //u: posicao da corda
	/*    u_old em t-1    */
	/*    u_now em t      */
	/*    u_new em t+1    */
    //k = dt/dx (aqui, dx = 1, k=dt)
    double u_new[jmax], t, u_old[jmax], u_now[jmax], k, erro, tmax;

	//tmax: tempo final
    tmax = 100;

	//variamos o dt para fazer um grafico do erro em função de dt
    int n=1;
    while( n < 1200)
    {
        k = 0.001*n;

		//condicao inicial
        for (j = 0 ; j < jmax ; j++)
        {
            u_now[j] = sin(M_PI*j/(jmax - 1));
        }
        
		//condicao de contorno
		u_new[0] = 0;
		u_new[jmax-1] = 0;
		
		u_old[0] = 0;
		u_old[jmax-1] = 0;
		
		//calculo de u para t=-dt, utilizando o metodo de leapfrog - aqui usamos que du/dt = 0 em t=0
        for (j = 1; j< jmax-1 ; j++)
        {
            u_old[j] = u_now[j] + 0.5 * pow(k,2) * (u_now[j+1] - 2 * u_now[j] + u_now[j-1]);
        }

		//atualizacao da onda
        for(t = 0 ; t < tmax ; t+=k)
        {
        	//para que o ultimo tempo em que se calcula a posicao seja tmax, independentemente do valor de dt
            if(t+k >tmax)
            {
                k = tmax - t;
            }
            
            atualizar_onda(u_new,u_now,u_old,jmax,k);

            memcpy(u_old,u_now, sizeof(double)*jmax);
            memcpy(u_now,u_new, sizeof(double)*jmax);
        }

		//calculo do erro
        erro = 0;
        for(j = 0 ; j < jmax ; j++)
        {
            double analitica = cos((M_PI*t)/(jmax - 1))*sin((M_PI*j)/(jmax - 1));
            erro += fabs(analitica - u_now[j]);
        }
        erro = erro/jmax;
        fprintf(arq,"%.12lf %.12lf\n", n*0.001,erro);


        n++;

    }

    fclose(arq);

}


void atualizar_onda(double u_new_[],double u_now_[],double u_old_[],int jmax,double k)
{
    int j;


    for(j = 1 ; j < jmax-1 ; j++)
    {
        u_new_[j] =   2 * u_now_[j] - u_old_[j] + pow(k,2) *( u_now_[j+1] - 2 * u_now_[j] + u_now_[j-1] );
    }


}