%% Matlab script 
%% Created by JKV, April 4. 1999
clear;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This script solves an equation of the type
%% -(a(x)*u'(x))' = f(x)  on [0,1]
%% with boundary conditions of one of the following types.
%% It uses the simplest type of FEM for an arbitrary division of [0,1].
%% The actual work is done in the accompanying function 'femsolver.m',
%% which expects to get the correct input arguments, as explained below.
%% Here the exact solution is expected to be in file 'u.m', but 
%% is only needed for plotting and error analysis.
%% Functions a(x) > 0, f(x) and x*f(x) should be given in function
%% .m-files which are passed to femsolver in the strings 
%% afu, ffu and xffu.
afu='aa'; ffu='ff'; xffu='xff';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Mixed natural boundary conditions ("Robin" or "Neumann"):
%% If you want to use boundary conditions of the type
%% -u'(0) + a0*u(0) = b0 and/or
%%  u'(1) + a1*u(1) = b1
%% choose lbctype = 1 and/or rbctype = 1 and specify the 
%% respective parameters below:
a0 = 0; b0 = 0;
a1 = 0; b1 = 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Essential boundary conditions ("Dirichlet"): 
%% If you want to use boundary conditions of the type
%% u(0) = u0 and/or 
%% u(1) = u1 
%% choose lbctype = 0 and/or rbctype = 0 and specify the
%% the values for u0 and/or u1
u0 = 0;
u1 = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Select boundary condition types. Set lbctype and rbctype to
%% the values 0 or 1 (independently of each other) to 
%% choose the boundary condition types to be used on the left
%% and right, respectively.
lbctype = 0; 
rbctype = 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Boundary condition info given above is passed to the femsolver 
%% function as vectors lbcparam and rbcparam. 
lbcparam(1)=lbctype; lbcparam(2)=a0; lbcparam(3)=b0; lbcparam(4)=u0;
rbcparam(1)=rbctype; rbcparam(2)=a1; rbcparam(3)=b1; rbcparam(4)=u1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Choose the node points 0=x(1),x(2),...,x(m+1),x(m+2)=1.    
%% These could be obtained as input from a file or typed from 
%% the keyboard or, as is done here, generated automatically. 
m=10;
x=(0:(1/(m+1)):1)';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Get the FEM-solution and plot it with the exact solution u(x).
sol=femsolver(x,lbcparam,rbcparam,afu,ffu,xffu);
plot(x,sol,'bd-'); hold on;
plot(x,u(x),'r-'); hold off;


%%%%%%%%%%%%%%%%%%%%
%% Error analysis %%
%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Here a uniform h=1/(m+1) is assumed on each round.

format long;
n=4;

for i=1:n,
  clear x sol exact 
  m(i)=5*2^(i-1)-1;
  h(i)=1/(m(i)+1);
  x=0:h(i):1;
  exact=u(x);
  sol=femsolver(x,lbcparam,rbcparam,afu,ffu,xffu);
  eh(i)=max(abs(sol-exact'));
end

ratio=eh(2:n)./eh(1:(n-1));

[h' eh' [1 ratio]']
