function x0new = TRBDFDiffusion(D,InVec,p,dt,nsteps)
%----------------------------------------------------------------------------------- 
% Time marching the diffusion step using the TR-BDF2 rule (stable for large dt)
% @ Sohan Dharmaraja MIT JUN2007


m=p^2;
dx = 1/(p+1);
DiffCoeff = D*dt/(4*(dx^2));
MyCoeff   = -(D*dt)/(2*(dx^2));
% f=10*ones(m,1)/(p+1)^2;

Ab=zeros(m,p+1); Ab(:,1)=4; Ab(:,2)=-1; Ab(p:p:m,2)=0; Ab(:,p+1)=-1;

% Solve using MATLAB
A=spdiags(Ab(:,[1,2,p+1]),-[0,1,p],m,m); A=A+tril(A,-1)';
x0 = InVec;
I = sparse(1:m,1:m,1);
A = DiffCoeff*A;
MMatrix = (I + A)\(I - A);

for i = 1:nsteps
    x0tmp = MMatrix*x0;
    x0new = (3*I - MyCoeff*A)\((4*I + MyCoeff*A)*x0tmp - x0);
    x0 = x0new;
end