#!/usr/bin/octave
disp("Tests for  u_t-u_{xx}=0 u(0)=u(pi)=0 u(0,x)=noncontinuous function ...
(characteristic function of [1,2]");
global L h mesh
#L approx of -laplacian*h^2 h space mesh step f 
function y=F(u,t=0)
global L h mesh

y=-(1/(h*h))*L*u + ff(mesh,t);
end

function y=ff(x,t=0)
 y=0*x;
end


M=100;
t=linspace(0,.01,30);
#diffusion
 h=pi/M;
 mesh=x=(1:M-1)'*h;
 L=(2)*sparse(diag(ones(M-1,1)))- sparse(diag(ones(M-2,1),1)) - sparse(diag(ones(M-2,1),-1));
 u0=(x<2)&(x>1);
 X=lsode(@F,u0,t);
 #plot(x,u0,x,X');pause(1);
 for j=1:length(t),
  title=[";diffusion t=",num2str(t(j)),";"];
  plot(x,u0,x,X(j,:)',title);pause(0.5);
 end
printf("press any key for longer time nad nonzero  force term\n"); 
pause;

#We can add force function f(t)=f equal one for x>2.5;  otherwise zero


function y=ff(x,t)
#f(x) = 1 x>2.5
#f(x) = 0  wpp
 y=(x>2.5);
end


t=linspace(0,2,60);
 grid on
 h=pi/M;
 x=(1:M-1)'*h;
 L=(2)*sparse(diag(ones(M-1,1)))- sparse(diag(ones(M-2,1),1)) - sparse(diag(ones(M-2,1),-1));
 u0=(x<2)&(x>1);
 X=lsode(@F,u0,t);
 #plot(x,u0,x,X');pause(1);
 for j=1:length(t),
  title=[";diffusion + force t=",num2str(t(j)),";"];
  plot([0; x; pi],[0; u0; 0],[0; x; pi],[0 X(j,:) 0]',title);pause(0.5);
 end

#compare with solution of -u''=f u(0)=u(pi)=0 - stationary state
disp("Error between stationary state and u(t,x) for t=2 ");
uinf=h*h*(L\ff(x,0));
norm(uinf- X(length(t),:)','inf')

pause(1);

disp("..and t = 100");
#for longer time
t=linspace(0,100);
X=lsode(@F,u0,t);
norm(uinf- X(length(t),:)','inf')

#error is like 1e-9 - we can assume that they are almost the same

