clear all
mu=0.01; tau=100; L=1; M=100; N=2000; x=linspace(0,L,M)';
v=sqrt(tau/mu); dx=x(2)-x(1); dt=dx/v; p=(v*dt/dx)^2;
%s=-2*x.^2 + 0.8*x + 0.92; g=transpose(sin((1:M)*2*3.14)); l=sin((0:N)'*dt*1000); r=-0.28*ones(length(1),N);
h=0.8 * ones(length(1),N);

s = exp(-5*(x-0.8).^2);
g = x.^2-x;
l = 0.04*ones(1,N);
k = -1.6*ones(1,N);

b = 0.5;
beta = b/2/mu;
q = 1+beta*dt;
u = 1-beta*dt;

f(:,1)=s;
%f(M,2)=p*(k(1)*2*dx + f(M-1,1) - 2*f(M,1) +f(M-1,1)) + +2/q*f(M,1)-u/q*(f(M,2) - g()*2*dt); % 13
f(M,2)=(2*p*(k(1)*dx  - f(M,1) +f(M-1,1)) + +2/q*f(M,1) + u/q*g(M)*2*dt)/(1+u/q); % 13
f(2:M-1,2)=p/2*(f(3:M,1)-2*f(2:M-1,1)+ f(1:M-2,1))+f(2:M-1,1)+u*dt*g(2:M-1);
f(1,2)=l(2);
for n=2:N-1
    f(1,n+1)=l(n+1);
    f(M,n+1)=2*p/q*(k(1)*2*dx + f(M-1,n) -2*f(M,n)- +f(M-1,n) )+2/q*f(M,n)-u/q*f(M,n-1); %13,15,17,18
    f(2:M-1,n+1)=p/q*(f(3:M,n)-2*f(2:M-1,n)+f(1:M-2,n))+2/q*f(2:M-1,n)-u/q*f(2:M-1,n-1);

     clf;
     plot(1:M,f(:,n));
     ylim([-4,4]);
     drawnow;
end


% f(:,1)=s;
% f(1,2)=p*(f(2,1)- f(1,1)-dx*h(1)) + f(1,1) +dt*g(1);
% f(2:M-1,2)=p/2*(f(3:M,1)-2*f(2:M-1,1)+ f(1:M-2,1))+f(2:M-1,1)+dt*g(2:M-1);
% f(M,2)=r(2);
% for n=2:N-1
%     f(1,n+1)=2*p*(f(2,n)-f(1,n)- dx*h(n) )+2*f(1,n)-f(1,n-1);
%     f(2:M-1,n+1)=p*(f(3:M,n)-2*f(2:M-1,n)+f(1:M-2,n))+2*f(2:M-1,n)-f(2:M-1,n-1);
%     f(M,n+1)=r(n+1);
%     clf;
%     plot(1:M,f(:,n));
%     ylim([-4,4]);
%     drawnow;
% end

surf(f)
shading interp

% plot(x, f(:,N),'LineWidth', 2);
%
% hold on
% yE = 0.8.*x - 1.08;
% plot(x, yE,'r','LineWidth', 2);
% grid on
%
% legend('last iteration', 'equilibrium');