% Initialization.
close all, clear all

l = 10; N = 101; dx = l/(N-1);
x0 = 5; v = 0.8; t0 = 5;
alpha = 0.1; eta = 0; gamma = -0.2;

A1 = -2.*eye(N) + diag(ones(N-1,1),1) + diag(ones(N-1,1),-1);
A1(1,2) = 2;
A1(N,N-1) = 2;
A = A1./dx^2;

Alpha = alpha.*eye(N);

Gamma = gamma.*ones(N,1);
Gamma(1,1) = gamma - 2*eta/dx;
Gamma(N,1) = gamma + 2*eta/dx;

% Initial values.
x = 0:dx:l;
c = (x-x0)/sqrt(1-v^2);
u0(1:N,1) = 4.*atan(exp(c));
u0(N+1:2*N,1) = -2*v/sqrt(1-v^2).*sech(c);

% Solution through ode45.
options = odeset('RelTol', 1e-8, 'AbsTol', 1e-8);
[tout,uout] = ode45(@sineG, [0 t0], u0, options, A, Alpha, Gamma, N);

% Exact solution.
uexact = 4.*atan(exp(c));

% Plot the solution.

%tsize = size(tout), ysize = size(yout), xsize = size(x)
%surf(x,tout,yout)

figure(1)
plot(x,uout(1:N),'-b',x,uexact,'-r')
xlabel('x'), ylabel('u(x,0)')
legend('Num. sol.', 'Exact sol.')