dam break simulation using lax wendroff and 1D st venant equation

i am trying to spproximate a PDE in the form below using the lax wendroff 2 step method in matlab coding:

[h ; hu ] = [ hu ; hu^2 + 1/2gh^2] = [0; -ghbz] (where bz will equal zero)

i believe this is then the case

d(h)/dt + d(hu)/dx = 0

and

d(hu)/dt + d(hu^2 + 1/2gh^2)/dx = 0

as the st venant equation was formed by depth deriving the 1D continuity equation

however i have attempted to program this and i keep produced a column matrix full of "NaN" could some one try to point out where i have been going wrong? thanks kyle

code for matlab:

%%%%%%%%%%%%%%%% function compare close all clc clear all

%intial values ntime = 1000000; dt=0.00050; nx=100; time = 0; a=2;output=0.24; %step size calculation dx= (1/nx); %create size of u_int vector u_int = zeros(nx,2); %create little u_int vector for the initial values of height begining and %end depending which value of A is used u_int(nx,2) =1; u_int(1,1) = 1;

%loop for the two directions needed for vec = 1:2 % %lax_wendroff [H] = lax_wendroff_2(nx,a,output, dt, ntime, time,u_int, vec,dx); figure(vec) H plot(H(:,vec),'r*');grid on;legend('centered'); hold on end end

function [H] = lax_wendroff_2(nx,a,output, dt, ntime, time,u_int, vec,dx); %compute original size matracies g = 9.81; uh=zeros(nx,2); update = uh ; H = uh; Hy = uh; UH = uh; if vec == 1; %for i = 1:nx u_int(nx-2:nx,2) = 1; u_int(1:3,1) = 1; %end for i = 1:nx UH(i,vec)=u_int(i,vec); Hy(i,vec)=u_int(i,vec); uh(i,vec)=u_int(i,vec); H(i,vec)=u_int(i,vec); A=a; end else uh(nx,vec)=u_int(nx,vec); h(nx,vec)=u_int(nx,vec); A=-a; end for p = 1:ntime; if(time + dt>output); dt=output-time; end % calculate coefficent coeff=dt/dx; %calculate interim steps for i = 1:(nx-1) % first derivative Hy(i+1,vec)= ((H(i+1,vec)+H(i,vec))/2) - ((coeff*A/2)*((uh(i+1,vec) - uh(i,vec))));

                    UH(i+1,vec)=((uh(i,vec)+uh(i+1,vec))/2) ;
                                -coeff*A*(((uh(i+1,vec)^2)/(H(i+1,vec)+g/2*H(i+1,vec)^2)) - ((uh(i,vec)^2)/(H(i,vec)+g/2*H(i,vec)^2)));

                end
                %calculate full time steps
                for i = 2:(nx-1)
                    update_h(i+1,vec) = H(i,vec) - coeff*A*(UH(i+1,vec) - UH(i,vec));%%%%check i notation and convention may be wrong way

                    update_uh(i+1,vec) = (uh(i,vec)) - (coeff*A*(((UH(i+1,vec)^2)/(Hy(i+1,vec)+g/2*Hy(i+1,vec)^2))- ((UH(i,vec)^2)/(Hy(i,vec)+g/2*Hy(i,vec)^2))));
                end
        for i = 3:(nx-1)
            H(i,vec) = update_h(i,vec);
            uh(i,vec) = update_uh(i,vec);
        end
            time = time + dt;
if time==output
break   
end

end end %%%%%%%%%%%%%%%%%