function [u, A, b] = MyPoissonSolver(p, t, e, EqData, BdryData)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% Purpose:  Solves D(-aDu)+cu=f in a domain D with
%%           boundary conditions anDu=k(g-u) where nDu
%%           nDu denotes the derivative out of D.
%%   
%%        p  1 x n vector of node coordinates.
%%
%%        t  3 x m matrix of element node numbers t(1,el)=nn1
%%           and t(2,el)=nn2, 1 <= el <= m, pointing to 
%%           corresponding node coordinates p(nn1) and p(nn2),
%%           and subdomain reference tags t(3,el), useful
%%           if data are given by different expresssions in
%%           different parts of the domain.
%%
%%        e  2 x q matrix of boundary node numbers and 
%%           reference tags, with node number e(1,bel)=nn pointing
%%           to the boundary node coordinate p(nn), and e(2,bel)
%%           is a reference tag (see example below).
%%
%%   EqData  matlab function returning [a,c,f] given
%%           the coordinate x and subdomain reference tag.
%%
%% BdryData  matlab function returning [k,g] given x and
%%           a boundary reference tag.
%% 
%%        A  n x n global stiffness matrix
%%
%%        b  n x 1 global load vector
%%
%%        u  n x 1 solution at node points p
%%           
%%
%% Example:  To solve D(-(1+x)Du)=1/x in [0,1] with boundary
%%           conditions Du(0)=u(0) and u(1)=7 (k = oo), define
%%           p = 0:0.01:1, t = [1:100; 2:101; ones(1, 100)],
%%           e = [1 101; 1 2] and functions
%%
%%           function [a, c, f] = EqData(x, tag)
%%             a = 1+x; c = 0; f = 1/x;
%%
%%           and
%%
%%           function [k, g] = BdryData(x, tag)
%%             if tag==1  k = 1; g = 0; end
%%             if tag==2  k = 10000; g = 7; end
%%
%%           To call the solver, type
%%           [u, A, b] = MyPoissonSolver(p, t, e, 'EqData', 'BdryData')
%%
%%           To plot the solution, type  plot(p, u)
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% initializing

  n = size(p, 2);
  A = zeros(n, n);
  b = zeros(n, 1);

% collecting element contributions to A and b

  for el = 1:size(t, 2)
    [dA, db] = ElementContributions(p, t, EqData, el);
    nn = t(1:2, el);
    A(nn, nn) = A(nn, nn) + dA;
    b(nn)     = b(nn)     + db;
  end

% collecting boundary contributions to A and b

  for bel = 1:size(e, 2)
     [dA, db] = BoundaryContributions(p, e, BdryData, bel);
     nn = e(1, bel);
     A(nn, nn) = A(nn, nn) + dA;
     b(nn)     = b(nn)     + db;
  end

% solving resulting system

  u = A\b;





% subroutines


function [dA, db] = ElementContributions(p, t, EqData, el)

% returns the element stiffness matrix dA and element load
% vector db for element number el defined by t(1:3,el) and
% p(t(1:2,el), and the equation data returned by EqData,
% using element midpoint quadrature.

  NodeCoords = p(t(1:2, el));
  x = sum(NodeCoords)/2;             % element midpoint!
  [a, c, f] = EqData(x, t(3, el));

  dx =
  phi =
  Dphi =

  dA =                % the 2x2 element stiffness matrix
                      % int_el a Dphi(I) Dphi(J) + c phi(I) phi(J) dx
                      % I,J = 1,2

  db =                % the 2x1 element load vector
                      % int_el f phi(J) dx, J=1,2



function [dA, db] = BoundaryContributions(p, e, BdryData, bel)

% returns the contributions dA and db from boundary element/node
% bel to the stiffness matrix and load vector.

  NodeCoord = p(e(1, bel));
  x = NodeCoord;
  [k, g] = BdryData(x, e(2, bel));
  ds = 1;
  phi = 1;

  dA = k*phi*phi*ds;     % the 1x1 bdry element contribution
                         % int_bel k phi(I) phi(J) ds, I,J = 1,2

  db = k*g*phi*ds;       % the 1x1 bdry element contribution
                         % int_bel k g phi(J) ds, J=1
  


function [a, c, f] = EqData(x, tag)

% returns the values at x of the data a=a(x), c=c(x) and f=f(x) in the
% differential equation D(-aDu)+cu=f

  a = 1;
  c = 0;
  f = 1/x;



function [k, g] = BdryData(x, tag)

% returns the values at x of the data k=k(x) and g=g(x)

  if tag==1
     k = 10000;
     g = 0;
  elseif tag==2
     k = 1;
     g = 1;
  end