nowait;
border(1,0,1,10) begin
      x:=t; y:=0; end;
border(1,0,2,20) begin
      x:=1; y:=t; end;
border(1,0,1,10) begin
      x:=1-t; y:=2; end;
border(1,0,2,30) begin
      x:=0; y:=2-t; end;
border(2,0,2*pi,30) begin
     x:= 0.4 + 0.2 * cos(-t);
     y:=0.5+ 0.2 * sin(-t);
end;
buildmesh(1000);

euler=0; crankn=0; gir=0;
mudt:= 0.05*0.1; i:=1; j:=2; k:=3;

iter(10) begin

  solve(euler,i) begin
      onbdy(1) euler=0;
      onbdy(2) euler=100;
      pde(euler) id(euler) - laplace(euler) 
                        * mudt  = euler;
  end;
 plot(euler); i:=-1;

  solve(u,j) begin
      onbdy(1) u=0;
      onbdy(2) u=100;
      pde(u) id(u)*2 - laplace(u) * mudt 
                       = crankn*2;
  end;
  crankn = 2*u-crankn;  plot(crankn); j:=-2;

  if(k==3) then oldgir = euler;
  solve(gir,k) begin
      onbdy(1) gir=0;
      onbdy(2) gir=100;
      pde(gir) id(gir)*3 - laplace(gir) * mudt 
                      = 2*gir + oldgir;
      oldgir=gir;/*assignment before solve*/
  end;
plot(gir); k:=-3;

end; /* end loop*/

save('gir.dta',gir);
save('crankn.dta',crankn);
save('euler.dta',euler);