function [x,v]=mypde3(N,M)

% same thing with CN 

T = 4; 
kappa = 1; 

k = T/M ;
h = 1/N ;
r = k * kappa/h^2 ; 
x = (0:N)'*h; 
t = (0:M)*k; 

u = zeros(N+1,M+1); 
u(:,1) = 0 ;             
u(1,:) = 0 ;         
u(N+1,:) = 2*sin(5*t) ;           

plot(x,u(:,1))
axis([0,1,-2,2])
drawnow 

a = (1+r)*ones(N-1,1);
b = (-r/2)*ones(N-2,1);
c = (-r/2)*ones(N-2,1);

for j=1:M
   R = (1-r)*u(2:N,j) + r/2*u(1:(N-1),j) + r/2*u(3:(N+1),j);
   R(1) = R(1) + r/2*u(1,j+1) ;
   R(N-1) = R(N-1) + r/2*u(N+1,j+1); 
   u(2:N,j+1) = tridiag(a,b,c,R);
   plot(x,u(:,j+1))
   axis([0,1,-2,2])
   drawnow
end

v = u(:,M+1);