% Az egységnégyzeten oldja meg a program a 
% u'_t-Laplace u =-20 hõvezetési egyenletet
% y=0 esetén u=1-x, x=0 esetén u'_x=y^2-1, 
% különben homogén Dirichlet-peremfeltétellel. 
% A kezdeti függvény (1-x)(1-y).

u0=@(x,y)(1-x).*(1-y);
f=@(x,y)-20+0*(x+y);

% A merevségi és tömegmátrix konstrukciója
n=15;
dt=0.01;
h=1/(n+1);
A=kron(gallery('tridiag',n+1),speye(n))+kron(speye(n+1),gallery('tridiag',n));
A(1:n,1:n)=A(1:n,1:n)/2;

D1=gallery('tridiag',n+1,0,0,1);
D2=gallery('tridiag',n,0,1,1);
D=kron(D1,D2);

B=(kron(speye(n+1),gallery('tridiag',n,1,6,1))+D+D')*h^2/12;
B(1:n,1:n)=B(1:n,1:n)/2;

% A terhelési vektor
[xi,yi]=meshgrid(0:h:1-h,1-h:-h:h);
b=h^2*reshape(f(xi,yi),n^2+n,1);
b(1:n)=b(1:n)/2;
b(1:n)=b(1:n)+h*(1-[1-h:-h:h].^2)';
b(2*n:n:n+n^2)=b(2*n:n:n+n^2)+(1-[h:h:1-h])';
b(n)=b(n)+1/2;

% A kezdeti függvény ábrázolása 

u=reshape(u0(xi,yi),n^2+n,1);

z=zeros(n+2);
z(2:n+1,1:n+1)=reshape(u,n,n+1);
z(n+2,1:n+2)=1-[0:h:1];
tri=[0 0 0];
for i=1:(n+1)*(n+2)
    if mod(i,n+2)~=0
        tri(length(tri(:,1))+1,:)=[i,i+1,i+n+3];
        tri(length(tri(:,1))+1,:)=[i,i+n+2,i+n+3];
    end
end
tri(1,:)=[];
[xi,yi]=meshgrid(0:h:1,1:-h:0);
trisurf(tri,xi,yi,z)
axis equal
axis([0,1,0,1,-2,1])


pause

% Idõlépések és a megoldás ábrázolása.

for k=1:1000
u=(B+dt*A)\(B*u+dt*b);

z=zeros(n+2);
z(2:n+1,1:n+1)=reshape(u,n,n+1);
z(n+2,1:n+2)=1-[0:h:1];
tri=[0 0 0];
for i=1:(n+1)*(n+2)
    if mod(i,n+2)~=0
        tri(length(tri(:,1))+1,:)=[i,i+1,i+n+3];
        tri(length(tri(:,1))+1,:)=[i,i+n+2,i+n+3];
    end
end
tri(1,:)=[];
[xi,yi]=meshgrid(0:h:1,1:-h:0);
trisurf(tri,xi,yi,z)
axis equal
axis([0,1,0,1,-2,1])
pause(0.00001)
end