import pymatbridge
ip = get_ipython()
pymatbridge.load_ipython_extension(ip)

%%matlab
path(path,'/Users/rjl/chebfun/chebfun_v4.2.2889/chebfun')

%%matlab
N = 100;
x = chebpts(N);
x = x(2:N);
Dop = chebop(@(u) diff(u));
D = Dop(N);
D = -D(2:N,2:N);



%%matlab
[V,Lam] = eig(D);
Lam = diag(Lam);
Lmax = max(abs(Lam))

%%matlab
nmax = 2000;  % number of time steps to reach t=1
dt = 1./nmax;
plot(real(Lam*dt), imag(Lam*dt),'r.')
eitheta = exp(1i*linspace(0,2*pi,1000));
hold on
plot(eitheta-1,'k')
axis 'equal'
xylim = dt*Lmax + 0.1;
axis([-xylim,0.1,-xylim,xylim])

%%matlab
x1 = linspace(-xylim,0.1,100);
y1 = linspace(-xylim,xylim,200);
[xx,yy] = meshgrid(x1,y1);
zz = xx + 1i*yy;
I = eye(size(D,1));
sigmin = zeros(size(zz));
for j=1:length(x1)
    for i=1:length(y1)
        sigmin(i,j) = min(svd(zz(i,j)*I - dt*D));
        end
    end

plot(real(Lam*dt), imag(Lam*dt),'r.')
hold on
plot(eitheta-1,'k')
axis 'equal'
axis([-xylim,0.1,-xylim,xylim])
contour(xx,yy,sigmin,[1e-2,1e-3,1e-4],'b')
contour(xx,yy,sigmin,[1e-5,1e-6,1e-7],'k')
contour(xx,yy,sigmin,[1e-8,1e-9,1e-10],'r')


%%matlab
u = exp(-50*(x+0.5).^2);
plot(x,u,'b');

%%matlab
for n=1:nmax
    u = u + dt*D*u;
    end

plot(x,u,'r-o')
hold on
utrue = exp(-50*(x-0.5).^2);
plot(x,utrue,'b')

%%matlab
norm(u-utrue,inf)

%%matlab
un = [0; u];
pn = chebfun(u)
phat = chebfun('exp(-50*(x-0.5).^2)')
max_error = norm(pn-phat)

%%matlab
un = exp(-50*(x+0.5).^2);
plot(x,un,'b');

dt = 1e-5;
nmax = 3000;
unm1 = exp(-50*(x+0.5-dt).^2);

for n=1:nmax
    unp1 = unm1 + 2*dt*D*un;
    unm1 = un;
    un = unp1;
    end
    
hold on
plot(x,un,'r-o')


%%matlab
k = 10;
plot(x,real(V(:,k)))
Lam(k)