function [x,U] = implicitpde(N,T)

   dx = (1-0)/N;
   x = dx*[0:N];
   dt = dx*dx/2.0;
   
   U(:,1) = cos(2*pi*x');
   A = diag((1+2*dt/(dx*dx))*ones(1,N-1))+diag(-dt/(dx*dx)*ones(1,N-2),1)+...
       diag(-dt/(dx*dx)*ones(1,N-2),-1);
   for n = 1:T,
      b = U(2:N,1);
      b(1) = b(1)+dt/(dx*dx);
      b(end) = b(end)+dt/(dx*dx);
      U(2:N,1) = A\b;
      U(1,1) = 1;
      U(N+1,1) = 1;
   end