% This code uses matlab's ode45 to numerically solve the wave
% equation.

% 2004-03-01 kfp

% integrate from xf to -xf 
xf = 50;

% initial condition (x +ve large) for each integration
% k^+ = 1/2, w = 1
kp=1/2;
km=1;

u0=[exp(i*kp*xf);i*kp*exp(i*kp*xf)];

% solve the ode specify in wave1_f.m
opt = odeset('AbsTol', atol);
opt = odeset(opt, 'RelTol', rtol);
[x, u] = ode45(@wave1_f, [xf -xf], u0, opt);

M=[exp(-i*km*xf), exp(i*km*xf);i*km*exp(-i*km*xf), -i*km*exp(i*km*xf)];
b=[u(end,1);u(end,2)];
z=M\b

%plot the wavespeed profile with a = 1,2,10
figure(1); clf; hold on;
y=-6:0.05:6;
c_1=(3+tanh(y))/2;
c_2=(3+tanh(2*y))/2;
c_3=(3+tanh(20*y))/2;
plot(y, c_1, 'b--');
plot(y, c_2, 'r-');
plot(y, c_3, 'k--');
xlabel('x');
ylabel('c(x)');

%phase plot of \eta(x) for both real and imaginery parts
figure(2);  clf;  hold on;
axis([-2 2 -2 2]);
plot(real(u0(1)), real(u0(2)), 'ro');
plot(real(u(:,1)), real(u(:,2)), 'b-');
xlabel('real \eta(x)-axis');
ylabel('real \eta''(x)-axis');   % note '' will render as '
title('real phase portrait');

figure(3);  clf;  hold on;
axis([-4 4 -4 4]);
plot(imag(u0(1)), imag(u0(2)), 'ro');
plot(imag(u(:,1)), imag(u(:,2)), 'b-');
xlabel('imaginery \eta(x)-axis');
ylabel('imaginery \eta''(x)-axis');   % note '' will render as '
title('imaginery phase portrait');

%plot eta(x) vs x
figure(4);  clf;  hold on;
plot(x, real(u(:,1)), 'b.-');
plot(x, imag(u(:,1)), 'r--');
xlabel('x');
ylabel('\eta');

T = 1/z(1)
R = z(2)*T
Texact = 2/3
Rexact = -1/3
Terr = abs(T-Texact)
Rerr = abs(R-Rexact)
