function [u]=ue(t,x,y)
    u=exp(-t)*cos(x+y);
endfunction
function [f]=fe(t,x,y,mu)
    f=-exp(-t)*cos(x+y)-mu*exp(-t)*(-cos(x+y)-cos(x+y));
endfunction


A=0.; B=0.5; C=0.; D=0.5; T=0.5; mu=2.5;

nx=16; ny=nx; dx=(B-A)/nx; dy=(D-C)/ny; dt=min(dx,dy);

for i = 1:nx+1
        for j = 1:ny+1
        u0(i,j)=ue(0.,A+(i-1)*dx,C+(j-1)*dy);
    end;
end; 

u1=u0; 

// time loop
n=0; t=0.; kmax=0;
while ((abs(T-t) > 10^(-12)) & (n<500)),
    //chose next time step
    dt=min(dt,dx,dy,T-t);

    //boundary condition plus initial guess for gauss-seidel iteration
    for i = 1:nx+1
        for j = 1:ny+1
        u1(i,j)=ue(t+dt,A+(i-1)*dx,C+(j-1)*dy);
        end;
    end;
    u1(2:nx,2:ny)=u0(2:nx,2:ny);
    
    c=1 +mu*dt/dx^2 +mu*dt/dy^2
    
    //gauss-seidel iterations
    dif=1.; k=0; ua=u1;
    while ((dif>2*dx^4) & (k<2500)),
        // interior points
        for i = 2:nx
        xi=A+(i-1)*dx;
        for j = 2:ny
            yj=C+(j-1)*dy;
            
            b=u0(i,j)+0.5*dt*(fe(t,xi,yj,mu)+fe(t+dt,xi,yj,mu));
            b=b+0.5*dt*mu*(u0(i-1,j)-2.*u0(i,j)+u0(i+1,j))/dx^2;
            b=b+0.5*dt*mu*(u0(i,j-1)-2.*u0(i,j)+u0(i,j+1))/dy^2;
           
            u1(i,j)=b/c;
            u1(i,j)=u1(i,j) +(0.5*dt*mu*(u1(i-1,j)+u1(i+1,j))/dx^2)/c;
            u1(i,j)=u1(i,j) +(0.5*dt*mu*(u1(i,j-1)+u1(i,j+1))/dy^2)/c;         
        end;
        end;        
               
        dif=max(abs(u1-ua)); k=k+1; ua=u1;
    end;
    if (k>kmax) then kmax=k; end;
    
    dif,k
        
    n=n+1; t=t+dt; u0=u1;
end
n,t 

for i = 1:nx+1
    for j = 1:ny+1
        u0(i,j)=ue(t,A+(i-1)*dx,C+(j-1)*dy);
    end;
end; 

printf('case= %d\n',nx);
printf('error= %.15f\n',max(abs(u1-u0)));
printf('max number of iterations: %d\n',kmax);
