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

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

%%matlab
beta = 50;
t0 = -1./(4*beta);
uhat = @(x,t) sqrt(-t0/(t-t0)) * exp(-x.^2 / (4*(t-t0)));
x = linspace(-3,3,1000);
t = 0.1;
plot(x,uhat(x,t),'b')
hold on
plot(x,0*x,'k')
plot([-1,-1],[-.25,.25],'k')
plot([1,1],[-.25,.25],'k')

%%matlab
ubar = @(x,t) uhat(x,t)-uhat(x-2,t)-uhat(x+2,t);
plot(x,ubar(x,t),'b')
hold on
plot(x,0*x,'k')
plot([-1,-1],[-.25,.25],'k')
plot([1,1],[-.25,.25],'k')

%%matlab
uhat(1,0.02)

%%matlab
x = chebfun('x');
L = chebop(@(u) diff(u,2));

%%matlab
Nvals = 20:5:40;  % spatial grids resolutions
for n=Nvals
    Ln = L(n);
    Ln = Ln(2:n-1,2:n-1);
    x = chebpts(n);
    x = x(2:n-1);
    errors = [];
    dts = [];
    
    % number of time steps logarithmically spaced:
    nmax_vals = floor(logspace(2,3.5,8));
    for nmax = nmax_vals
        tfinal = 0.02;
        dt = tfinal/nmax;
        dts = [dts, dt];
        A = eye(n-2) - 0.5*dt*Ln;
        
        % initial data:
        u = exp(-beta*x.^2);
        
        % time stepping loop:
        for nt=1:nmax
            u = A\(u + 0.5*dt*Ln*u);
            end

        utrue = ubar(x,tfinal);
        error = norm(u-utrue,inf);
        errors = [errors; error];
        end

    loglog(dts,errors,'o-')
    hold on
    text(3e-6,errors(end),sprintf('n = %g',n),'fontsize',15)
    end

% plot triangle showing slope dt^2 for expected second order accuracy 
xs1 = 2e-5;
xs2 = 1e-4;
ys1 = 1e-8;
ys2 = (xs2/xs1)^2 * ys1;
xs = [xs1,xs2,xs2,xs1];
ys = [ys1,ys1,ys2,ys1];
loglog(xs,ys,'r')
set(gca,'fontsize',15)
title('Log-log plot of error vs. \Delta t')