
clear

n = 99;
X = linspace(0,1,n+2);
h = X(2)-X(1);
k = 0.01;
c = 1;

U0 = max(0,min(10*X-5,-10*X+7));
v0 = zeros(n+2,1);
U = U0';
U(:,2) = U(:,1)+k*v0+c^2*k^2/2*zeros(n+2,1);

A = -1/h^2*(2*diag(ones(n,1))-diag(ones(n-1,1),-1)-diag(ones(n-1,1),1));

T = 0:k:1.2;

[T,V] = ode15s(@(t,v) [zeros(n,n) eye(n,n);A zeros(n,n)]*v,T,[U0(2:end-1)';zeros(n,1)]);
U = [zeros(1,length(T));V(:,1:n)';zeros(1,length(T))];
W = [zeros(1,length(T));V(:,n+1:end)';zeros(1,length(T))];

figure(1)
mesh(X,T,U')
axis tight
view(20,45)
xlabel('$x$','interpreter','latex')
ylabel('$t$','interpreter','latex')
zlabel('$u(x,t)$','interpreter','latex')
title('Solution $u(x,t)$','interpreter','latex')
set(gca,'fontsize',15)

figure(2)
mesh(X,T,W')
axis tight
view(20,45)
xlabel('$x$','interpreter','latex')
ylabel('$t$','interpreter','latex')
zlabel('$u_t(x,t)$','interpreter','latex')
title('Speed $u_t(x,t)$','interpreter','latex')
set(gca,'fontsize',15)

figure(3)
T = 0:k:5;
[T,V] = ode15s(@(t,v) [zeros(n,n) eye(n,n);A zeros(n,n)]*v,T,[U0(2:end-1)';zeros(n,1)]);
U = [zeros(1,length(T));V(:,1:n)';zeros(1,length(T))];
W = [zeros(1,length(T));V(:,n+1:end)';zeros(1,length(T))];

for i=1:length(T)
    plot(X,U(:,i),'linewidth',1)
    drawnow
    set(gca,'xlim',X([1 end]),'ylim',[-1 1]);
    M(i) = getframe;
end
