

n=10   
h=1/n #rom
m=100
T=2
k=T/m #tid

r=k/h^2

x=0:h:1;



%u=zeros(n+1,n+1,m+1);
u(2:end-1,2:end-1,1)=1;

%eksplisitt
for j=1:m
%   u(end,:,j)=sin(pi*x);
   for i=2:n
       for k=2:n
           u(i,k,j+1)=u(i,k,j)+.1*r*(u(i+1,k,j)+u(i,k+1,j)-4*u(i,k,j)+u(i-1,k,j)+u(i,k-1,j));
       end
   end
end


%implisitt
% A=(1+2*r)*diag(ones(n-1,1))-r*diag(ones(n-2,1),-1)-r*diag(ones(n-2,1),1);
% A=kron(eye(n-1),A);
% A=A-r*diag(ones((n-1)^2-n+1,1),n-1)-r*diag(ones((n-1)^2-n+1,1),n-1);
% u=ones((n-1)^2,m+1);


% for j=1:m 
%     u(:,j+1)=A\u(:,j);
% end


%crank-nicolson
% r=r/2;
% A_v=(1+2*r)*diag(ones(n-1,1))-r*diag(ones(n-2,1),-1)-r*diag(ones(n-2,1),1);
% A_h=(1-2*r)*diag(ones(n-1,1))+r*diag(ones(n-2,1),-1)+r*diag(ones(n-2,1),1);
% 
% for j=1:m 
%        u(2:n,j+1)=A_v\(A_h*u(2:n,j));
% end

figure


for j=1:m
    surf(x,x,u(:,:,j))
    axis([0 1 0 1 -1 1.2])
    pause(.1)
end
