#!/usr/bin/octave
disp("Tests for  u_t-u_{xx}=0 u(a)=u(b)=0 u(0,x)=sin(x) ; u(t,x)=exp(-t)sin(x)");
disp("FDM discretizations on equidistant mesh");
function y=F(x,t)
#vector field function for lsode
#user befora callin F(x,t) has to define global variables:
#L tridiagonal symmetric matrix (M-1)*(M+1) - 2 on diagonal -1 on sub- superdiagonals 
#h=(b-a)/M - mesh parameter (space discretization)
#x - vector of size M-1
global L h 
y=-(1/(h*h))*L*x;
end



global  L h
t=[0,1,2,4,5,100];
em=e0=zeros(length(t)-1,5);
M0=M=10;
%sin(x) on [0,1]
a=0;b=pi;U0=@sin; %function handle to initial value function
for k=1:5,
 hh(k)=h=(b-a)/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=U0(x);
 X=lsode(@F,u0,t);
 for j=2:length(t),
   em(j-1,k)=norm(X(j,:)'-exp(-t(j))*u0,'inf'); #max norm for t(j)
   e0(j-1,k)=sqrt(h)*norm(X(j,:)'-exp(-t(j))*u0,2);
 endfor  
 M*=2;
endfor
e=em;#discrete max norm
e=e0;#discrete L2 norm
function []=printerrors(em,e0,t,hh,M0,type=0)
M=M0;
#printing results
if(type==0)
 e=em;#discrete max norm
 printf("\n\nErrors in max discrete norm\n");
else
 e=e0;#discrete L2 norm
 printf("\n\nErrors in L2 discrete norm\n");
endif
printf("h:                ");
for k=1:5,
  printf(" %6.3e",hh(k));
 M*=2;
endfor
printf("\n");
for j=1:(length(t)-1),
 printf("t=%4.3e error ",t(j+1));
 for k=1:5,
  printf(" %6.3e",e(j,k));
 endfor
 printf("\n            ratio    ");
 for k=1:5,
  if(k>1)
     printf("%6.3f    ",e(j,k-1)/e(j,k));
     else printf("       ");
  endif 
 endfor
  printf("\n");
endfor
end %end printerror()
printerrors(em,e0,t,hh,M0,type=0);
printerrors(em,e0,t,hh,M0,type=1);
disp("t=100 - too long time interval");

%solving u_t=u_{xx} on [0,3] with zero bnd and 
%inital value a characteristic function of (1,2) 
a=0;b=3;U0=@(x) (x<2)&(x>1);
M=60;
h=(b-a)/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=U0(x);
 t=linspace(0,2,10);
 X=lsode(@F,u0,t);
for k=1:length(t),
 plot(x,X(k,:)); pause(1);
end
mesh(X);pause(1);

%same but for u0(x)=0 except u0(1)=10 or 100 etc
a=0;b=3;
M=60;
h=(b-a)/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=zeros(size(x));
 u0(20)=100;
 t=linspace(0,2,10);
 X=lsode(@F,u0,t);
for k=1:length(t),
 plot(x,X(k,:)); pause(1);
end
mesh(X);

