% Crank Nicolson solver for u_t = kappa u_{xx} + F(u)  0<x<L, 0<t<T  
%        with Dirichlet boundary conditions 
% using full Newton to solve at each step 

% parameters of the problem 
T = 8;
L = 4;
kappa = 0.1;

% set up lattice
M = 3000;
N = 100;
h = L/N;
k = T/M;
r = k*kappa/h^2;
x = (0:N)'*h;
t = (0:M)*k;
u = zeros(N+1,M+1); 

% boundary conditions and initial condition
u(1,:) = 0 ;
u(N+1,:) = 0 ;
u(:,1)  = 0;
u(N/2,1) = 1; 
plot(x,u(:,1))
axis([0,L,0,1])
drawnow 
pause

% numerical stuff 
a0 = (1+r)*ones(N-1,1);     % basic tridiagonal matrix
b  = -r/2*ones(N-2,1); 
c  = -r/2*ones(N-2,1); 

for j=1:M 

    U = u(:,j); 
    V = u(:,j+1);
    FU = F(U); 
    V(2:N) = U(2:N) + r*(U(3:(N+1))-2*U(2:N)+U(1:(N-1))) + k*FU(2:N);    % Euler to start
    FV = F(V); 

    e = V(2:N)-U(2:N) - r/2*(U(3:(N+1))-2*U(2:N)+U(1:(N-1))+...
               V(3:(N+1))-2*V(2:N)+V(1:(N-1))) - k/2*(FU(2:N) + FV(2:N)) ; 
    a = a0 - k/2*Fp(V(2:N)); 

    while norm(e)>1e-12
        delta = tridiag(a,b,c,e);
        V(2:N) = V(2:N) - delta ;
        FV = F(V); 
        e = V(2:N)-U(2:N)-r/2*(U(3:(N+1))-2*U(2:N)+U(1:(N-1))+V(3:(N+1))-2*V(2:N)+V(1:(N-1)))-k/2*(FU(2:N) + FV(2:N)) ; 
        a = a0 - k/2*Fp(V(2:N));
    end
    
    u(:,j+1)=V; 

    plot(x,u(:,j+1))
    axis([0,L,0,1])
    drawnow 
end