#!/usr/bin/octave
printf("\nDirichlet bnd condition\n\n");
#testing -u''+cu=f u(a)=al u(b)=be

f=@(x) sin(x);a=al=be=c=0;b=pi;N=10;sol=@sin;
#solution u(x)=sin(x)
errm=errl2=0;
printf("\n-u''=sin u(0)=0=u(pi) solution  sin(x)\n");
for k=1:6,
  errmp=errm;errl2p=errl2;
  [u,x]=FDsolver(f,c,a,b,al,be,N);
  errm=norm(u-sol(x),'inf');
  errl2=norm(u-sol(x),2)/sqrt(N);
  N*=2;
  printf("N=%d max %7.3f  %5.5f  L2 %7.3f  %5.5f \n",N,errm,errmp/errm,errl2,errl2p/errl2);
end 
#Convergence order  - two
printf("Convergence order - two\n");

#same for -u''=u u(0)=u(b)=1 b=20

f=@(x) 0;c=1; a=0;b=20; al=be=1;N=10; sol=@(x) (exp(x-20)+exp(-x))/(exp(-20)+1);
errm=errl2=0;
printf("\n-u''=u u(0)=1=u(20) solution  (exp(x-20)+exp(-x))/(exp(-20)+1)\n");
for k=1:6,
  errmp=errm;errl2p=errl2;
  [u,x]=FDsolver(f,c,a,b,al,be,N);
  errm=norm(u-sol(x),'inf');
  errl2=norm(u-sol(x),2)/sqrt(N);
  N*=2;
  printf("N=%d max %7.3f  %5.5f  L2 %7.3f  %5.5f \n",N,errm,errmp/errm,errl2,errl2p/errl2);
end 
#Convergence order  - two
printf("Convergence order - two\n");

printf("\nMixed bnd condition\n\n");
#testujemy -u''+cu=f u'(a)=al u(b)=be
printf("\n -u''=-2 u'(0)=0=u(1) solution  x*x\n");
f=@(x) -2;a=c=0;al=0;b=1;be=b*b;N=10;sol=@(x) x.*x;
#solution u(x)=x*x
errm=errl2=0;
for k=1:6,
  errmp=errm;errl2p=errl2;
  [u,x]=FDmixlft(f,c,a,b,al,be,N);
  errm=norm(u-sol(x),'inf');
  errl2=norm(u-sol(x),2)/sqrt(N);
  N*=2;
  printf("N=%d max %7.3f  %5.5f  L2 %7.3f  %5.5f \n",N,errm,errmp/errm,errl2,errl2p/errl2);
end 
#Convergence order - one
printf("Convergence order - one\n");

printf("\n-u''=sin u'(0)=1 u(pi)=0 solution  sin(x)\n");
f=@(x) sin(x);a=c=0;al=cos(a);b=pi;be=0;N=10;sol=@(x) sin(x);
#solution u(x)=sin(x)
errm=errl2=0;
for k=1:6,
  errmp=errm;errl2p=errl2;
  [u,x]=FDmixlft(f,c,a,b,al,be,N);
  errm=norm(u-sol(x),'inf');
  errl2=norm(u-sol(x),2)/sqrt(N);
  N*=2;
  printf("N=%d max %7.3f  %5.5f  L2 %7.3f  %5.5f \n",N,errm,errmp/errm,errl2,errl2p/errl2);
end 
#Lepiej - Convergence order  - two? Why?
printf("Convergence order - two!? Why? \n");


#same for -u''=u u'(0)=(exp(-20)-exp(-0))/(exp(-20)+1)   u(b)=1 b=20

f=@(x) 0;c=1; a=0;b=20; al=(exp(-20)-exp(-0))/(exp(-20)+1) ;
be=1;N=10; sol=@(x) (exp(x-20)+exp(-x))/(exp(-20)+1);
errm=errl2=0;
printf("\n-u''=u u'(0)=1=u(20) solution  (exp(x-20)+exp(-x))/(exp(-20)+1)\n");
for k=1:6,
  errmp=errm;errl2p=errl2;
  [u,x]=FDmixlft(f,c,a,b,al,be,N);
  errm=norm(u-sol(x),'inf');
  errl2=norm(u-sol(x),2)/sqrt(N);
  N*=2;
  printf("N=%d max %7.3f  %5.5f  L2 %7.3f  %5.5f \n",N,errm,errmp/errm,errl2,errl2p/errl2);
end 
#Convergence order  - one
printf("Convergence order - one\n");
