				%initial/boundary conditions
				%Lax-Wendroff method for u_t + au_x = 0

h = input('Please enter a value for h');
rho = input('Please enter a value for rho');
k = rho*h;
a = 1;

%setting max time, max space values
xmax = 5;
tmax = 5;

%obtaining x and t values from 0 to xmax and 0 to tmax respectively
x = 0:h:xmax;
t = 0:k:tmax;

%this values gives a benchmark for which space step results in x=1
%to be used for inputting the initial condition
step_ic = 1/h+1;

%Max space steps and time steps 
N = length(x)-1;
J = length(t);

%need extra inital information in order to use Lax-Wendroff 
%for non-periodic boundary conditions.  This give the base
%of domain over which computations actually occur
xextra = xmax + h*(J+1);
xlong = 0:h:xextra;
Nextra = length(xlong)-1



  %use lax-wendroff to compute u over trapezoidal domain (domain
  %of numerical dependence needed for given tmax

u = zeros(Nextra+1,J+1);

%initial condition
u(1:1:step_ic,1) = x(1:1:step_ic);
u(step_ic:1:2*step_ic-1,1) = 2-x(step_ic:1:2*step_ic-1);
u(2*step_ic-1:1:Nextra+1,1) = 0;


for j=1:1:J
  %setting left boundary values to zero
u(1,j+1) = 0;

u(2:1:Nextra-j+1,j+1) = (1-rho^2*a^2).*u(2:1:Nextra-j+1,j)+.5*rho*a*(rho*a-1).*u(3:1:Nextra-j+2,j) + .5*rho*a*(rho*a + 1).*u(1:1:Nextra-j,j);
end


%truncate to keep only solution u in domain where we want solution,
% i.e., 0<=x<=xmax, 0<=t<=tmax;
u_mod = zeros(N+1,J+1);
u_mod(:,:) = u(1:1:N+1,:);


%plotting solutions at t=0,1,2,3,4,5
figure;
for tstep = 1:(1/k):J+1
time = (tstep-1)*k
tstep
subplot(2,3,time+1)
plot(x,u_mod(:,round(tstep)))
end


%find true_u solution to compare to computed u_mod and obtain error results
true_u = zeros(N+1,6);


for n = 0:1:5

  i = round(n/k + 1)
  delta_x = round(2/h + 1) 
  %  u(x,t) = u_0(x-t)
    if i + delta_x -1 <= N+1
     true_u(i:1:i+delta_x-1,i) = u(1:1:delta_x,1);
    else 
      true_u(i:1:N+1,i) = u(1:1:N+2-i,1);
    end

       error(:,n+1) = u_mod(:,i)- true_u(:,i);
end

%print out error at times t=0,1,2,3,4,5
error
