function v=difusao7(t,x)
%function v=difusao7(t,x)
% objetivo eh resolver problema de difusao
%     v_t = nu * v_xx
%     v(x,t) = f(x)
%     v(0,t) = a(t), v(1,t) = b(t) , t > 0
% usando o esquema de Crank-Nicholson
% entrada:
%   t = [t0 t1 ... tN] : discretizacao no tempo
%   x = [x0 x1 ... xK] : discretizacao no espaco
% saida:
%   v : matriz com a solucao discretizada
% funcoes internas:  f(x), a(t) , b(t)
%  Mai 2012   
nu=10^(-3);

dt=t(2)-t(1); N=max(size(t))-1;

dx=x(2)-x(1); K=max(size(x))-1;

v = zeros(N+1,K+1);  % inicializa a matriz resposta

for k=1:K+1
  xk = (k-1)*dx; v(1,k)=f( xk ); % condicoes iniciais
end
for i=1:N+1
  ti =  (i-1)*dt; v(i,1) = a( ti ); v(i,K+1)=b(ti); % C.C.
end
r = dt*nu/(dx*dx); rs2=r/2;
vw=zeros(K-1,1);   % vetor do lado direito (b)
A = zeros(K-1,K-1); % matriz de coeficientes (A)
A(1,1)=1+r;
for j=2:K-1
  A(j,j-1) = -rs2; A(j,j) = 1+r; A(j-1,j) = -rs2;
end
% construimos a matriz A
R = chol(A); % notacao do Matlab: R'*R = A
% assim, para resolver A u = w :
% R'*R* u = w <==>  z = (R')\w ; u = R \ z
for i=1:N
  for k=2:K
     vw(k-1) = v(i,k) + rs2*( v(i,k+1)-2*v(i,k)+v(i,k-1) );
  end
  vw(1) = vw(1) + rs2*v(i+1,1); % condicao de contorno a esq
  vw(K-1) = vw(K-1) + rs2*v(i+1,K+1); % C.C a direita
  vz = (R')\vw; vu = R \ vz;  % A vu = vw
 
  v(i+1,2:K)=vu'; % salva a resposta na estrutura de dados
end


function u=f(x)
u = 1+3*sin(2*pi*x);
return
function u=a(t)
 u = 1;
return
function u=b(t)
 u = 1;
return


return
