Erro: Lax-Wendroff de dois passos

#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_laxw.dat", "w");

    int tmax,i,j,jmax;
    
    //tamanho da corda: jmax-1
    jmax = 50;
    
    //u: posição 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: tempo final
    tmax = 100;


	//variamos o dt para fazer um gráfico do erro em função de dt
    int n=1;
    while( n < 1200)
    {
        k = 0.001*n;
		
		//condição inicial
        for (j = 0 ; j < jmax ; j++)
        {
            u_now[j] = sin(M_PI*j/(jmax - 1));
        }

		//condição de contorno
		u_old[0] = 0;
		u_old[jmax-1] = 0;
		
		u_new[0] = 0;
		u_new[jmax-1] = 0;
		
		//cálculo de u para t=-dt, utilizando o método 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]);

        }


		//atualização da onda
        for(t = 0 ; t < tmax ; t+=k)
        {
        	//para que o último tempo em que se calcula a posição 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);
        }

		//cálculo 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;

	//para j=1 e j=jmax-2 (penúltimas posições), necessitamos de u em x=-1 e x=jmax, respectivamente, sendo que essas posições ficam fora do vetor (vai de 0 a jmax-1). Aqui, foi feito u(x=-1)=-u(x=1) e u(x=jmax)=-u(x=jmax-2)
	j=1;
	
	u_new_[j] = 2 * u_now_[j] - u_old_[j] + 0.5 * pow(k,2) * (0.5 * u_old_[j+2] - u_old_[j+1] + u_old_[j] - u_old_[j-1] - 0.5 * u_old_[j] + u_now_[j+1] - 2 * u_now_[j] + u_now_[j-1]);

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



}