function v = difusao2d_1(t,x,y)
//
//  x=[0 0.2 0.4 0.6 0.8 1]; y=[0 0.25 0.50 0.75 1.0]; t=[0 0.2 0.4 0.6];
nu=1/6;
dt = t(2)-t(1); N = max(size(t,1),size(t,2))-1;
dx = x(2)-x(1); J = max(size(x,1),size(x,2))-1;
dy = y(2)-y(1); K = max(size(y,1),size(y,2))-1;

rx = nu*dt/(dx*dx); ry=nu*dt/(dy*dy); // parametros

v=zeros(J+1,K+1,N+1); // array 3-dimensional

v(:,:,1)=feval(x,y,f); // condicao inicial

for i=1:N
   //
   // fronteiras em x=0 e x=1:
   for k=1:K+1
     v(1,  k,i+1) = g(x(1),y(k),t(i+1));
     v(J+1,k,i+1) = g(x(J+1),y(k),t(i+1));
   end
   // fronteiras em y=0 e y=1:
   for j=1:J+1
     v(j,  1,i+1) = g(x(j),  y(1),t(i+1));
     v(j,K+1,i+1) = g(x(j),y(K+1),t(i+1));
   end
   // agora no interior, usando o esquema numerico
   // implementacao mais simples: varredura bi-dimensional
   for j=2:J
     for k=2:K
       delta2_x = v(j+1,k,i)-2*v(j,k,i)+v(j-1,k,i);
       delta2_y = v(j,k+1,i)-2*v(j,k,i)+v(j,k-1,i);

       v(j,k,i+1)=v(j,k,i) + rx*delta2_x + ry*delta2_y + dt*F(x(j),y(k),t(i));
     end
   end
end  
// use mesh(y,x,v(:,:,i)) para visualizar
endfunction

//  abaixo: as funcoes da condicao inicial e condicoes de contorno
//  cuidado: nao esqueca que devem satisfazer EQ de COMPATIBILIDADE !!

function u=f(x,y)
  u = 1 + 2*x*(1-x)*y*(1-y)
endfunction

function u=g(x,y,t)
  u = 1 + sin( x*(1-x) ) * sin( y*(1-y) );
endfunction

function u=F(x,y,t)
  u = 0;
endfunction
 
