clear all; close all;

% Shock

for (m = 1:4)
    n = 8*m;
    h = 4/n;
    dt = h/2;
    
    u(1,1:4*n+1,m) = [ones(1,2*n) 0.5 zeros(1,2*n)];
    for (i = 2:n+1)
            for (j = 2:4*n)
                u(i,j,m) = ( (u(i-1,j-1,m))^2 - (u(i-1,j+1,m))^2 )/4/h*dt + (u(i-1,j+1,m)+u(i-1,j-1,m))/2;
            end
    end

    for (k = 1:2*n+1)
        y(k,m) = u(n+1,k+n,m);
    end 
end

figure(1);
x = linspace(-4, 4, 16+1);
subplot(4,1,1); plot(x, y(1:length(x),1)); legend('gird = 8'); xlabel('x'); ylabel('u'); title('Shock');
x = linspace(-4, 4, 16*2+1);
subplot(4,1,2); plot(x, y(1:length(x),2)); legend('gird = 16'); xlabel('x'); ylabel('u');
x = linspace(-4, 4, 16*3+1);
subplot(4,1,3); plot(x, y(1:length(x),3)); legend('gird = 24'); xlabel('x'); ylabel('u');
x = linspace(-4, 4, 16*4+1);
subplot(4,1,4); plot(x, y(1:length(x),4)); legend('gird = 32'); xlabel('x'); ylabel('u');


for (k = 1:n+1)
    t(k) = 2*dt*k;
end

v = zeros(n+1,2*n+1);
for (i = 1:n+1)
    for (j = 1:2*n+1)
        v(i,j) = u(i,j+n,4);
    end
end

figure(2);
title('Shock');
mesh(x,t,v);
xlabel('x'); ylabel('t'); zlabel('u');


%Fan
for (m = 1:4)
    n = 8*m;
    h = 4/n;
    dt = h/2;
    
    u(1,1:4*n+1,m) = [zeros(1,2*n) 0.5 ones(1,2*n)];
    for (i = 2:n+1)
            for (j = 2:4*n)
                u(i,j,m) = ( (u(i-1,j-1,m))^2 - (u(i-1,j+1,m))^2 )/4/h*dt + (u(i-1,j+1,m)+u(i-1,j-1,m))/2;
            end
    end

    for (k = 1:2*n+1)
        y(k,m) = u(n+1,k+n,m);
    end 
end

figure(3);
x = linspace(-4, 4, 16+1);
subplot(4,1,1); plot(x, y(1:length(x),1)); legend('gird = 8'); xlabel('x'); ylabel('u'); title('Fan');
x = linspace(-4, 4, 16*2+1);
subplot(4,1,2); plot(x, y(1:length(x),2)); legend('gird = 16'); xlabel('x'); ylabel('u');
x = linspace(-4, 4, 16*3+1);
subplot(4,1,3); plot(x, y(1:length(x),3)); legend('gird = 24'); xlabel('x'); ylabel('u');
x = linspace(-4, 4, 16*4+1);
subplot(4,1,4); plot(x, y(1:length(x),4)); legend('gird = 32'); xlabel('x'); ylabel('u');

for (k = 1:n+1)
    t(k) = 2*dt*k;
end

v = zeros(n+1,2*n+1);
for (i = 1:n+1)
    for (j = 1:2*n+1)
        v(i,j) = u(i,j+n,4);
    end
end

figure(4);
title('Fan');
mesh(x,t,v);
xlabel('x'); ylabel('t'); zlabel('u');